Algorithm for Linear Programming

CHAPTER

2

A Primal-Dual Interior Point Algorithm for Linear Programming Masakazu Kojima, Shinji Mizuno, and Akiko Yoshise

Abstract. This chapter presents an algorithm that works simultaneously on primal and dual linear programming problems and generates a sequence of pairs of their interior feasible solutions. Along the sequence generated, the duality gap converges to zero at least linearly with a global convergence ratio (1 - Yf/n); each iteration reduces the duality gap by at least Yf/n. Here n denotes the size of the problems and Yf a positive number depending on initial interior feasible solutions of the problems. The algorithm is based on an application of the classical logarithmic barrier function method to primal and dual linear programs, which has recently been proposed and studied by Megiddo.

§1. Introduction This chapter deals with the pair of the standard form linear program and its dual. (P)

(D)

Minimize c T x subject to

XES

= {x E Rn : Ax = b, x

~

O}.

Maximize b T Y subject to (y,z)

E

T = {(y,z)

E

Rm+n: ATy

+Z=

c, Z

~O}.

Here R n denotes the n-dimensional Euclidean space, cERn and bERm are constant column vectors, and A is an m x n constant matrix. Throughout the chapter we impose the following assumptions on these problems: (a) The set Sf = {x E S : x> O} of the strictly positive feasible solutions of the primal problem (P) is nonempty. N. Megiddo (ed.), Progress in Mathematical Programming © Springer-Verlag New York Inc. 1989

2. Primal-Dual Interior Point Algorithm

30

(b) The set Tl = {(y,z) E T: Z > O} of the interior points of the feasible region T of the dual problem (D) is nonempty. (c) Rank A = m. By the well-known duality theorem, we know that both of the problems (P) and (D) have optimal solutions with a common value, say s*. Furthermore, it is easily seen that the sets of the optimal solutions of(P) and (D) are bounded. To solve the primal problem (P), we consider the problem of minimizing the objective function c T x with the additional logarithmic barrier term - J.L Ij=llog Xj over the relative interior SI of the feasible region S: Minimize c T x - J.L I" log Xj j=l

subject to x

E

SI.

Here J.L is a penalty parameter. The logarithmic barrier function method was used for linear programs by Frish [2J and recently by Gill et al. [4J to develop a projected Newton barrier method in connection with Karmarkar's projective rescaling method [8,9]. Megiddo [14J has made a deep analysis on the curve consisting of the solutions x(J.L) of (Pit) (J.L > 0). He has shown some nice duality relations on the logarithmic barrier function method for a symmetric pair of primal and dual linear programs. A new algorithm we propose here will be founded on his analysis and duality relations. But we will deal with the nonsymmetric pair of the primal and dual problems, (P) and (D). So we will briefly explain them in terms of the nonsymmetric pair. We first note that the objective function of (Pit) is convex. Hence each problem (Pit) has at most one local minimal solution. Furthermore, we can prove under assumptions (a) and (b) that each subproblem (Pit) has a unique global minimal solution x (see Section 2 of Megiddo [14J), which is completely characterized by the Karush-Kuhn-Tucker stationary condition (see, for example, Mangasarian [13J):

Ax - b

= 0,

X>O.

Here X = diag(xl, ... ,x,,) denotes the diagonal matrix with the diagonal elements Xl' ... , x" and e E R" the vector of l's. Introducing a slack variable vector Z = p.X- 1 e, we will write the stationary condition as

Zx - p.e Ax - b

ATy

=

0,

= 0,

+z -

c

X> 0,

(2.1)

= 0,

where Z = diag(z I ' ... , z"). The first equality can be further written component wise as Xi Z i

= P.

(i

=

1,2, ... ,n).

31

§l. Introduction

That is, all the products of the complementary pair of the variables Xi and Zj = 1,2, ... , n) take a common value Il, the value of the penalty parameter. Let (X(Il), Y(Il), z(ll» denote a unique solution of the system (2.1) for each Il > O. Then each point (Y(Il), Z(Il» lies in the interior of the dual feasible region T, that is, (Y(Il), z(Il» E rl. The gap between the primal objective value at X(Il) E SI and the dual objective value at (Y(Il), Z(Il» E r l turns out to be

(i

= (c T

c T X(Il) - b T Y(Il)

-

Y(Il)T A)X(Il)

= Z(Il)T X(Il)

= eTZ(Il)x(ll) = nil. System (2.1) also gives a necessary and sufficient condition (the KarushKuhn-Tucker stationary condition) for (Y(Il), Z(Il» to be a maximal solution of the subproblem: Maximize b T Y + Il subject to (y, z)

E

n

L log Zj

j=l

rl.

This subproblem can be regarded as an application of the barrier function method to the dual problem (D). See Theorem 3.1 of [14] for more detailed relations on the subproblems (P1l) and (DIl ). Let r denote the curve consisting of the solutions ofthe system (2.1), that is,

r

=

{(X(Il), Y(Il), z(Il»: Il >

OJ.

As we take the limit Il--+ 0, the curve r leads to a pair of optimal solutions x* of problem (P) and (y*,z*) of problem (D), which satisfies the strong complementarity such that zt = 0 if and only if xt > 0 (Proposition 3.2 of [14]). We have observed that if (X(Il), Y(Il), Z(Il» is a solution of the system (2.1) or if (x (Il), Y(Il), Z(Il» lies on the curve r then X(Il)

E

(2.2)

SI,

(Y(Il), Z(Il»

E

r

I ,

(2.3)

and (2.4)

Thus, tracing the curve r, which has been proposed in [14], serves as a theoretical model for a class of primal-dual interior point algorithms for linear programs. The new algorithm, which will be described in Section 2, can be regarded as a practical implementation ofthis model using Newton's method. Newton's method (or projected Newton's method) has been used in several interior point algorithms for linear programs (Gill et al. [4], Imai [6], Iri and Imai [7], Renegar [16], Yamashita [21]). In those algorithms it is applied to minimization of barrier functions or penalty functions, whereas in our algorithm it will be applied to the system of equations (2.1) that has been derived as the Karush-Kuhn-Tucker stationary condition for problem PIl , the

32


minimization of the logarithmic barrier function over the interior of the primal feasible region. Given an initial point (XO, yO, ZO) E SI X TI, the algorithm generates a sequence {(x k , y\ Zk)} in SI x TI. To measure a "deviation" from the curve r at each (Xk, y\ Zk) E SI X TI, we introduce the quantity (2.5)

where

f/ J.

= x~z~

k min --

(i

= 1,2, ... ,n),

k mIn {r J i •••I -- 1, 2, ••• , n },

(2.6)

•

(the average of hk'S). Obviously, we see that n k ~ 1 and that (Xk, y\ Zk) E SI X TI lies on the curve r if and only if n k = 1. When the deviation nO at the initial point (XO, yO, ZO) E SI X TI is large, we reduce not only the duality gap but also the deviation. Keeping a small deviation, we then direct our main effort to reducing the duality gap. Several properties of the global convergence of the algorithm will be investigated in Section 3. The algorithm has two control parameters. Specifically, we will see that if we make a suitable choice of the parameters then the generated sequence {(x\ y\ Zk) E SI X TI} satisfies the inequalities CTXk+l _ b Ty k+1

~

(1 _ 2/(nn k »(c Tx k

nk+l _ 2 ~ (1 - l/(n nk+l ~ 2

ifnk

+ l»(n k -

~

_

bTyk),

2) if n k > 2,

2.

(Corollary 2.2). The first inequality ensures that the duality gap decreases monotonically. By the second, we see that the deviation n k gets smaller than 3 within at most O(n log nO) iterations, and then the duality gap converges to zero linearly with the convergence rate at most (1 - 2/(3n». Each iteration requires O(n3) arithmetic operations. Generally an initial primal and dual feasible solution (XO, yO, ZO) E SI X TI is not known, so we have to incorporate a phase 1 procedure. This will be explained in Section 4, where we provide artificial primal and dual linear programs with a known feasible solution. System (2.1) gives some significant insights into the interior point algorithms developed so far (e.g., Adler et al. [1], Gay [3], Gill et al. [4], Imai [6], Iri and Imai [7], Karmarkar [8,9], Renegar [16], Todd and Burrell [19], Yamashita [21], Ye and Kojima [22J). Originally these algorithms were designed to apply to various different forms oflinear programs. But when they are adapted to the standard form linear program (P), most of them have the common property that the search direction along which new interior feasible solutions are generated from a current feasible solution x is composed of

§2. Algorithm

33

the two directions, X P X c and X Pe, although the coefficients associated with these two directions vary in different algorithms. Here P denotes the orthogonal projection matrix onto the subspace {w E R n : AXw = O}, X = diag(xl' x 2 , ••• , x n ) an n x n diagonal matrix, and e the n-dimensional vector of l's. This fact has recently been pointed out by Yamashita [21]. If (x,y,z) lies on the curve r consisting of the solutions of system (2.1), then these two directions coincide with the tangent vector of the curve r. This indicates that the local behaviors of these algorithms and ours are quite similar in neighborhoods of the curve r. We refer to the paper [15] by Megiddo and Shub, who have studied the boundary behavior of Karmarkar's projective rescaling algorithm and the linear rescaling algorithm. In the discussions above, the matrix X serves as a linear (or affine) rescaling that transforms the current feasible solution x to the vector e in the primal space. When we deal with the dual problem (D), the diagonal matrix Z-1 = diag(1/z1' 1/z 2 , ••• , l/zn) plays this role [1]. It is also interesting to note that the rescaling matrix X for the primal problem (P) has been utilized to generate feasible solutions of the dual problem (D) to get lower bounds for the primal objective function in some of the algorithms listed above [3, 19,22]. By the first equality of the system (2.1) we see that if (x, y, z) E r, then the rescalings by the matrix X in the primal space and by the matrix Z-1 in the dual space differ only by a constant multiple so that they are essentially the same. In our algorithm the common rescaling matrix (XZ- 1 )1/2, which is the geometric mean of the primal rescaling matrix X and the dual rescaling matrix z-1, will be applied simultaneously to the primal and dual spaces even if the point (x, y, z) does not lie on the curve r.

§2. Algorithm We begin by deriving the Newton direction of the system of equations (2.1) at (Xk, y\ Zk) E Sf X Tf. If we denote the right-hand side of equations (2.1) by H(x, y, z), then the Newton direction (~x, ~y, ~z) at (x\ y\ Zk) is defined by the system of linear equations H(Xk,Y\Zk) - DxH~x - DyH~y - DzH~z = 0, where DxH, DyH, and DzH denote the Jacobian matrices of the mapping H at (x\ y\ Zk) with respect to the variable vectors x, y, and z, respectively. It follows that Z~X

+ X ~z = v(fJ.), A~x

AT~y

where

= 0,

+ ~z = 0,

(2.7)

34


x = diag(x1, x~, ... ,x!), z = diag(zL z~, ... , z!), v(lt) = XZk - Ite E W.

It should be noted that a solution of this system depends linearly on the penalty parameter It, so we will write a solution as (~x(It), ~y(It), ~z(It)). By a simple calculation, we obtain ~x(lt) = ~Y(It)

{Z-l - Z-l XAT(AZ- 1XATr 1 AZ- 1 }v(It),

= -(AZ- 1 XAT)-l AZ-1v(It),

~z(lt) =

AT(AZ- 1XAT)-l AZ-1v(It).

We now introduce the diagonal matrix

D = (XZ- 1)1/2

= diag«x1/z1)1/2, ... , (x!/z!) 1/2), the vector-valued linear function (2.8)

or (i = 1,2, ... , n)

in the variable It, and the orthogonal projection matrix

Q = DAT(AD2 AT)-l AD onto the subspace {DA Ty:y can be written as

E

Rm}. Then the Newton direction (~x(It), ~Y(It),

~z(lt))

~x(lt) =

D(I - Q)u(It),

~Y(It) = -(AD2 ATr 1 ADu(It),

(2.9)

~z(lt) = D- 1 Qu(It)·

These equalities and (2.7) imply D- 1 ~x(lt)

o=

+ D~z(lt) = u(It),

~x(ltf ~z(lt) = (D- 1 ~x(It)f(D~z(It))·

(2.10)

At the start of the algorithm described below, we fix two real parameters (J and, such that 0 ~ , < (J < 1, which control the penalty parameter It and the step length, respectively. Define nk, f/ (i = 1,2, ... , n), fa~e, and f~in by (2.5) and (2.6). We will set the penalty parameter It = (Jfa~e and then compute a new point, the (k + 1)th point from (x k, yk, Zk) in the direction -(~x(It), ~Y(It), ~z(It)). That is, (Xk+1, yk+1, Zk+1) can be written in the form

35

§2. Algorithm

X(IX, Jl)

=

y(lX, Jl)

= yk

Xk - IXAx(Jl),

(2.11)

- IXAY(Jl),

Z(IX, Jl) = Zk - IXAz(Jl).

Here IX denotes a step length. As the step length IX changes, the product J; = XiZ i of each complementary pair of variables Xi and Zi changes quadratically: J;(IX, J1) = Xi (IX, Jl)Zi(lX, Jl)

= xfzf -

lX(xf AZi(Jl)

+ zfAxi(Jl)) + AXi(Jl)Azi(Jl)1X2

= J;k - (J;k - Jl)1X + AXi(Jl)Azi(J1)1X2

(i

= 1,2, ... , n), (2.12)

(since xfAzi(J1) + zfAxi(Jl) = J;k - Jl by (2.7)) and the average fave of J; = 1,2, ... , n) as well as the duality gap c T X - b T Y changes linearly:

(i

fave(lX, Jl) =

{~ J;(IX, J1) }/n

= x(a., Jlf Z(IX, Jl)/n = fa~e

- (fa~e - J1)1X

(since AX(Jlf AZ(Jl)

=

0 by (2.10)) (2.13)

CTX(IX, J1) - b Ty(lX, J1) = X(IX, Jlf Z(IX, Jl) =

XkT Zk - (XkT Zk - nJ1)1X

=

(CTX k - bTyk) - {(CTX k - bTyk) - nJl}IX. (2.14)

Figure 2.1 illustrates the functions J; and fave. If we neglected the quadratic term AX i(J1)Az i(J1)1X 2 in the function J;, it would induce a line through the two points (O,J;k) and (1, Jl). It should be noted that the latter point does not depend on the coordinate i. The step length IXk will be determined by f3k

= max {IX: J;(IX', Jl)

~ l/!(IX')

for all IX' E (0, IX) and i = 1,2, ... , n},

(2.15)

where (2.16) See Figure 2.2. Now we are ready to state the algorithm. Algorithm

Step 0: Let (XO,yO,ZO) E SI X TI be an initial point and G a tolerance for duality gap. Let (J and r be real numbers such that 0 ~ r < (J < 1. Let k = O. Step 1: c T Xk - b T yk < G then stop.

36


jj(a, Il)

fav.(a., Il)

J;(a., Il)

Figure 2.1. The function.t; (i = 1,2, ... , n) and f.ve.

Step 2: Step 3: Step 4: Step 5:

Define fa~e and f~in by (2.6). Let Ii = (Jfa~e and compute (~X(Ii), ~Y(Ii), ~Z(Ii)} by (2.9). Compute pk and cx k by (2.15). Let Xk+l

=

x(cxk, Ii),

yk+l

=

y(cx k, Ii),

Zk+l

=

z(cxk, Ii),

(2.17)

where x(cx, Ii), y(cx,,u), and z(cx,,u) are given by (2.11). Step 6: Let k = k + 1. Go to step 1.

§3. Global Convergence Properties In this section we show some global convergence properties of the algorithm. The following lemma provides a common lower bound for the functions /; (i = 1,2, ... , n) defined by (2.12).

37

§3. Global Convergence Properties

_-+_____

......L.._ _ _~-_+--~

l/!(rx)

J.(rx, (Jfa~c)

Figure 2.2. Computation of pk.

Lemma. Let Then h(a., /1-) PROOF.

~

((IX, /1-)

for any a. E (0, l)(i = 1,2, ... , n).

By (2.12) we have h(IX, /1-)

= hk

- (h k - /1-)a.

+ (D- I ~x(/1-)MD~z(/1-))ia.2.

The first two constant and linear terms on the right-hand side are bounded by f:'in - (f:'in - /1-)a. for all a. E [0,1]. So we have only to evaluate the last quadratic term. The vectors D- I ~x(/1-) and D~z(/1-) satisfy (2.10). In other words, the vector u(/1-) is the orthogonal sum of these two vectors. Hence (D- I ~x(/1-)MD~z(/1-))j

= (u(/1-)j -

(D~z(/1-n)(D~z(/1-))j

~ IU(/1-)jI2 /4

38


for alIj

=

1,2, ... , n. Therefore

(D- l AX(Jl»);(DAz(Jl)}i

Theorem 2.1 If the step length less than 1, then OCk

L (D-

= (D- l AX(JlW(DAz(Jl)) -

OCk

j¥oi

AX(Jl)MDAz(Jl»j

l

determined at step 4 of the kth iteration is

~ 4(u - .) k 2 ~ 4(u -2') k' - n(l - 2u + n u ) - n(l + u )n

CTXk+l _ bTI+l nk+l -

u/t

= (1 - (1 -

U)OCk)(CTXk -

~ (1 - v)(nk - u/.)

nk+l ~

u/.

if n k

(2.18)

bTyk),

(2.19)

nk,

(2.20)

if u/.
O"/r throughout all the iterations then O"/r < n k ~ (1 - v)k(nO - O"/r)

+ O"/r

(k

= 0,1, ... ).

In both cases we have that

nk

~

max{O"/r,nO}

(k = 0,1,2, ... )

and that n k gets smaller than (O"/r:) + 1 in at most O(n log nO) iterations, say iterations. Now assume that k ~ k. Then we have, by inequality (2.18), (

k

1 - O")a

k

(1 - 0")4(0" - r)

~ n(1 + 0"2)«0"/r) + 1)·

By inequality (2.19), the duality gap c TXk - b Tyk attains the given accuracy e and the iteration stops in at most O(n log«c TXO - b TyO)/e) additional iterations. We summarize the above results in the following corollary. Corollary 2.1. The algorithm terminates in at most O(n log nO)

+ O(n log«c TXO

- b TyO)/6))

iterations.

In practical implementations of the algorithm, there may be various ways of setting the control parameters 0" and r such that 0 ~ r < 0" < 1. Although the values of these parameters are fixed at the start of the iteration of the

41

§3. Global Convergence Properties

algorithm and the fixed values are used throughout the iterations, we could change those values at each iteration. As a special case, we obtain: Corollary 2.2. Suppose that

(J'

= 1/2 and. = 1/4. Then

IY. k ~ 4/(mr k), c TXk+1 _ b Tyk+1

~

(1 _ 2/(nnk)(c TXk _ b Tyk),

n k+ 1 - 2

~

(1 - 1/(1

nk+1

~

2 if n k ~ 2.

+ n))(n k -

2)

if 2 < nk,

PROOF. We have only to show the third inequality because the others follow immediately from Theorem 2.1. By the first inequality ofthe corollary, we see nklY. k ~ 4/n. Hence, substituting this inequality, (J' = 1/2, and • = 1/4 into (2.28), we obtain the desired result.

In view of the second inequality of Corollary 2.2, the kth iteration reduces the gap between the lower bound c TXk and the upper bound b Tyk for a common optimal value of (P) and (D) by at least 2/(nn k). The following theorem shows a method of getting better lower and upper bounds for the optimal value at each iteration; if we set Jl = 0 and take the largest step size IY. satisfying the primal and dual feasibility conditions x(IY.,O) E Sand (Y(IY.,O), Z(IY., 0)) E T, then we can reduce the gap by at least 1/(nnk)1/2. Theorem 2.2. Define

O=max{IY.:};(IY.',O)~O algorithm. Let

x=

for alllY.'E[O,IY.] (i= 1,2, ... ,n)} at step 3 of the

x(O,O),

y = y(O,O),

z = z(O,O).

Then x and (y, z) are feasible solutions of (P) and (D), respectively, and the inequality holds.

By the construction of the mappings x, y, and z whose definition has been given by (2.11), x and (.p, z) satisfy the equality constraints Ax = band ATY + z = c. Since };(IY., 0) = xi(lY., O)Zi(lY., 0) for every IY. ~ 0 (i = 1,2, ... , n) by definition (2.12), the nonnegativity x ~ 0 and z ~ 0 follows from the definition of O. Thus we have shown the first statement of the theorem. By inequality (2.14) with f.1. = 0 and IY. = 0, we have PROOF.

o ~ cTx -

bTy = (1 - O)(CTXk - bTyk).

(2.29)

This implies that 0 ~ 1. Applying the lemma as in the proof of Theorem 2.1, we then see that 0 is not less than a positive solution of the equation


42

It follows that

o ~ 2{ -f::'in + «(f::'in)2 + f::'inllu(O)112)1/2}/llu(O)112 =

2/{(1

+

lIu(O)1I2/f::'in)1/2

+ I}.

On the other hand, we see by definition (2.8) of u(ll) that lIu(O) 112

Hence

= nfa~e =

nf::'innk.

o ~ 2/{(1 + nnk)1/2 + I} ~

1/(nnk)1/2.

By substituting this inequality into (2.29), we obtain the desired inequality. Remark. If (x\ y\ Zk) lies on the curves r consisting of the solutions of the system (2.1) then n k = 1. In this case, the Newton direction (Ax(O), Ay(O), Az(O» defined by (2.9) with Il = 0 coincides with the tangent vector of the curve r. Theorem 2.2 above ensures that the linear extrapolation of the curve r along the tangent vector to the boundary of the feasible regions in the primal and dual spaces generates a pair of primal and dual feasible solutions that reduces the interval containing the common optimal value of the problems (P) and (D) by at least n- 1/ 2 •

§4. Artificial Primal and Dual Linear Programs for Initiating the Algorithm In the discussions so far, we have been assuming that a feasible solution E SO of the primal problem (P) and a feasible solution (yO, ZO) E TO of the dual problem (D) are available so that we can immediately start the algorithm. In general, however, it is necessary to provide a means for determining an (x, y, z) E SO x TO so that the algorithm can be initiated. The method described below is essentially due to Megiddo [14]. See Section 3 of [14]. Let (XO, yO, ZO) E R n + m + n be an arbitrary point satisfying XO > 0 and ZO > O. Let K and A be sufficiently large positive numbers, whose magnitude will be specified below. Corresponding to the pair of primal and dual linear programs, (P) and (D), which we want to solve, we consider the following pair of artificial primal and dual linear programs: XO

Minimize (Pi)

subject to

Ax (ATyO

+ ZO -

C)TX

+ (b -

AxO)x n +1

= b,

+ Xn + 2 =

A,

§4. Artificial Primal and Dual Linear Programs for Initiating the Algorithm

where

X n +1

and

X n +2

43

are artificial real variables.

Maximize subject to (D')

Ym+l

+ Zn+l = X, + Zn+2 = 0,

where Ym+l, Zn+l, and Zn+2 are artificial real variables. We need to take K and A. satisfying

A. > (ATyO

+ ZO - cfxo

(2.30)

and

(2.31) so that (XO, X~+1' X~+2) and (yO, y~+1' zo, Z~+l' Z~+2) are feasible solutions of(P') and (D'), respectively, where

X~+2

= A. - (AT yO + ZO

y~+1

= -1,

-

C)TxO,

Therefore, we can apply the algorithm to the artificial problems (P') and (D') with these initial feasible solutions. Let x* and (y*, z*) be optimal solutions of the original problems (P) and (D), respectively. The theorem below gives a sufficient condition on the constants K and A. for the algorithm applied to the artificial problems (P') and (D') to succeed in generating approximate solutions of the original problems (P) and (D). Theorem 2.3. In addition to (2.30) and (2.31), suppose that

A. > (AT yO

+ ZO -

C)T x*

(2.32)

and K

> (b - AxOfy*.

(2.33)

Then the following (a) and (b) hold:

x

(a) A feasible solution (x, Xn+1' Xn+2) of (P') is minimal if and only if is a minimal solution of (P) and Xn+1 = o. (b) A feasible solution CP,Ym+1,Z,Zn+1Zn+2) of (D') is maximal if and only if (y, z) is a maximal solution of (D) and Ym+1 = o.

44


PROOF.

Define

X:+2

= A. -

(ATyO

+ ZO -

C)TX*.

Then (x*, X:+l' x:+ 2) is a feasible solution. Let (x, x n +1, x n +2) be an arbitrary feasible solution of (PI) with X n +1 > 0. Then we see c T x*

+ KX:+ 1 =

y*Tb

+ (b - AxO)xn +1 } z*fx + I(X n + 1 (by Ay* + z* = C, Xn +1 > 0, and (2.33» + I(Xn +1 (by (Z*)T x ~ 0).

= yH {Ax

< (c ~

c TX

°

This implies that (x, x n +1 , x n +2) with Xn +1 > cannot be a minimal solution of (PI). In other words, any minimal solution of (PI) must satisfy X n +1 = 0. By the continuity, we also see that (x*, x:+1' x:+ 2) is a minimal solution of (PI). Therefore, if a feasible solution (x, xn +1 , Xn +2) of (PI) is minimal then xn +1 = 0 and c T X = c T x*. Since x satisfies all the constraints of (P), it must be a minimal solution of (P). Conversely, if (x, 0, xn + 2 ) is a feasible solution of (PI), and if x is a minimal solution of (P), then the objective value c T X + I(X n + 1 coincides with the minimal value c T x* + KX:+ 1. Hence it is a minimal solution of (PI). Therefore we have shown (a). We can prove (b) similarly. If we take XO = e, yO = 0, and ZO = e then the conditions (2.30), (2.31), (2.32), and (2.33) imposed on the constants I( and A. tum out to be A. > (e - C)T e = n - c T e,

(2.30)'

I( > 0,

(2.31 )'

A. > (e - C)TX*

and I(

= eTx* - cTx*,

> (b - Aefy* = bTy* - yH Ae.

(2.32), (2.33)'

§5. Concluding Remarks Kojima [10] has proposed a method of determining optimal basic variables in Karmarkar's projective rescaling algorithm. The idea of his method is based on the relaxation technique and the duality theorem (see Sections 5 and 6 of Kojima and Tone [12]), and could be incorporated in the algorithm to determine both primal and dual optimal basic variables. The system of equations (2.1) which has played a central role in the development of the algorithm can also be derived from the notion of a center (an analytic center) of a polytope, which has been proposed by Sonnevend [17] as a fundamental tool for a class of interior point algorithms for smooth

§5. Concluding Remarks

45

convex minimization over polyhedra including linear programs. Since the set of the optimal solutions of (P) is bounded by the assumptions (a) and (b) of the Introduction, the polyhedron S(s)

=

{(]l, x)

E

Rl+n:

X E

S,cTx +]l

= S,]l

~

O}

is bounded and contains a strictly positive (]l, x) for each s greater than the optimal value s* of (P). Hence, according to [17J, we can consistently define the center of S(s) to be a unique maximal solution of the problem Maximize

n

L: log + log]l Xj

j~l

subject to (]l, x) E S(s), (]l, x) > O. Using the Karush-Kuhn-Tucker optimality condition for this problem, we can easily verify that (]l, x) > 0 is a solution of (ps ) for some s > s* if and only if there exists a (y, z) E R m + n for which the system (2.1) is satisfied. Similarly, the center of the polytope T(t)

=

{(]l,y, z) E R m +n +1 : (y, z) E T, b T Y - ]l

=

t,]l ~

O},

which is induced from the dual problem (D), is defined to be a unique maximal solution of the problem Maximize

m

L: log(zi) + log]l

i~l

subject to (]l, y, z)

E

T(t), (]l, z) > O.

Then (]l, y, z) is a maximal solution of (D,) for some t < s* if and only if there exists an x E R n for which the system (2.1) is satisfied. Therefore the duality relations (Megiddo [14J) which we briefly explained in the Introduction could be stated in terms of the problems (Ps ) (s > s*) and (D,) (t < s*). Tracing the path of centers of the polytope T(t) has been employed in the polynomial time algorithm proposed by Renegar [16J and the multiplicative penalty function method by Iri and Imai [7J (see also Imai [6J). After writing the draft of this chapter, the authors learned that Tanabe [18J has developed an interior point algorithm based on the idea of tracing the path of centers in the primal and dual spaces. This work was presented at the conference on "Progress in Mathematical Programming" held in Monterey, California, on March 1-4, 1987. The authors learned there that, besides the papers referred in the Introduction and the comments above, many studies have been done on the paths of centers. Among others, Gonzaga [5J and Vaidya [20J both proposed algorithms that trace the path of cent~rs to solve linear programs in O(n 3 L) arithmetic operations. The algorithm given here requires O(n4 L) arithmetic operations. The algorithm will be modified and extended in [l1J so as to solve a class of linear complementarity problems including linear and convex quadratic programs in O(n3 L) arithmetic operations.

46


Acknowledgment The authors wish to thank Professor Masao Mori for his warm encouragement of their study on new interior point algorithms for linear programming and Dr. Nimrod Megiddo for inviting one of them to the conference on "Progress in Mathematical Programming" where this work was presented.

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17]

[18]

I. Adler, N. Karmarkar, M. G. C. Resende, and G. Veiga, An implementation ofKarmarkar's algorithm for linear programming, Working Paper, Operations Research Center, University of California, Berkeley (May 1986). K. R. Frish, The logarithmic potential method of convex programming, unpublished manuscript, University Institute of Economics, Oslo, Norway (1955). D. M. Gay, A variant of Karmarkar's linear programming algorithm for problems in standard form, Math. Programming 37 (1987),81-90. P. E. Gill, W. Murray, M. A. Saunders, J. A. Tomlin, and M. H. Wright, On projected Newton barrier methods for linear programming and an equivalence to Karmarkar's projective method, Math. Programming 36 (1986),183-209. C. C. Gonzaga, An algorithm for solving linear programming problems in O(n3 L) operations, in Progress in Mathematical Programming, N. Megiddo, (ed.), Springer-Verlag, New York, this volume. H. Imai, Extensions of the multiplicative penalty function method for linear programming, Journal of the Operations Research Society of Japan 30 (1987), 160-180. M. Iri and H. Imai, A multiplicative barrier function method for linear programming, Algorithmica 1 (1986),455-482. N. Karmarkar, A new polynomial-time algorithm for linear programming, Proc. 16th Annual ACM Symposium on Theory of Computing, Washington, D.C. (1984). N. Karmarkar, A new polynomial-time algorithm for linear programming, Combinatorica 4 (1984), 373-395. M. Kojima, Determining basic variables of optimal solutions in Karmarkar's new LP algorithm, Algorithmica 1 (1986),499-515. M. Kojima, S. Mizuno, and A. Yoshise, A polynomial-time algorithm for a class oflinear complementarity problems, Math. Programming (to appear). M. Kojima and K. Tone, An efficient implementation of Karmarkar's new LP algorithm, Research Report No. B-180, Department of Information Sciences, Tokyo Institute of Technology, Meguro, Tokyo (April 1986). O. L. Mangasarian, Nonlinear Programming, McGraw-Hill, New York, 1970. N. Megiddo, Pathways to the optimal set in linear programming, this volume. N. Megiddo and M. Shub, Boundary behavior of interior point algorithms in linear programming, to appear in Math. Oper. Res. (1988). J. Renegar, A polynomial-time algorithm, based on Newton's method, for linear programming, Math. Programming 40 (1988), 59-94. G. Sonnevend, An 'analytic center' for polyhedrons and new classes of global algorithms for linear (smooth, convex) programming, Proc. 12th IFIP Conference on System Modeling and Optimization, Budapest (1985), to appear in Lecture Notes in Control and Information Sciences, Springer-Verlag, New York. K. Tanabe, Complementarity-enforcing centered Newton method for linear programming: Global method, presented at the Symposium, or "New Method

§5. Concluding Remarks

[19] [20] [21] [22]

47

for Linear Programming," Institute of Statistical Mathematics, Tokyo (February 1987). M. J. Todd and B. P. Burrell, An extension of Karmarkar's algorithm for linear programming used dual variables, Algorithmica 1 (1986),409-424. P. M. Vaidya, An algorithm for linear programming which requires O«(m + n) x n2 + (m + n)1.5 n)L) arithmetic operations, AT&T Bell Laboratories, Murray Hill, N.J. (1987). H. Yamashita, A polynomially and quadratically convergent method for linear programming, Working Paper, Mathematical Systems Institute, Inc. Tokyo (October 1986). Y. Ye and M. Kojima, Recovering optimal dual solutions in Karmarkar's polynomial time algorithm for linear programming, Math. Programming 39 (1987),305-317.