- Email: [email protected]
An interior point method for solving systems of linear equations and inequalities
Operations Research Letters 27 (2000) 101–107
www.elsevier.com/locate/dsw
An interior point method for solving systems of linear equations and inequalities Markku Kallio ∗ , Seppo Salo Helsinki School of Economics, Runeberginkatu 14-16, Helsinki 00100, Finland Received 1 March 1999; received in revised form 1 March 2000
Abstract A simple interior point method is proposed for solving a system of linear equations subject to nonnegativity constraints. The direction of update is deÿned by projection of the current solution on a linear manifold deÿned by the equations. Infeasibility is discussed and extension for free and bounded variables is presented. As an application, we consider linear c programming problems and a comparison with a state-of-the-art primal–dual infeasible interior point code is presented. 2000 Elsevier Science B.V. All rights reserved. Keywords: Interior point algorithms; Linear inequalities; Linear programming
1. The problem and a method of solution An inspiration to our article arises from the vast literature on interior point methods for solving mathematical programming problems (see e.g. [2]). In particular, we are motivated by the success of primal–dual infeasible interior point algorithms for linear programming. Our aim is to solve systems of linear equations subject to nonnegativity constraints, for which optimality conditions for linear programming is a special case. We propose a simple interior point algorithm which can be initiated with any positive vector. It always converges, and if the problem is feasible, the limit point is a solution.
Let q ∈ Rm and Q ∈ Rm×n with m ¡ n and rank(Q) = m. Our problem is to ÿnd z ∈ Rn such that Qz = q; z¿0:
Corresponding author. E-mail address: kallio@hkkk.ÿ (M. Kallio).
(1.2) n
Denote M = {z ∈ R | Qz = q}; Z = diag(z), and for d ∈ Rn , deÿne Z −2 -norm by kdkZ −2 = (dT Z −2 d)1=2 . Our method starts with z 0 ∈ Rn ; z 0 ¿ 0. For iteration k; k = 0; 1; 2; : : : ; denote Z = diag(z k ), suppressing k in Z. The direction dk of update is deÿned by dk = z − z k where z is the point of M which is closest to z k in Z −2 -norm; i.e., z is the unique solution to the following problem: min kz − z k kZ −2 : z∈M
∗
(1.1)
(1.3)
A step in this direction is taken. If z k + dk ¿0 a solution for (1.1) – (1.2) is found; otherwise the step size k ¡ 1 is applied so that z k+1 = z k + k dk ¿ 0.
c 2000 Elsevier Science B.V. All rights reserved. 0167-6377/00/$ - see front matter PII: S 0 1 6 7 - 6 3 7 7 ( 0 0 ) 0 0 0 3 6 - 5
102
M. Kallio, S. Salo / Operations Research Letters 27 (2000) 101–107
The steps of the algorithm are formalized as follows: 0. Initialization. Set the iteration count k to 0, pick any z 0 ∈ Rn such that z 0 ¿ 0, and let the step size parameter  ∈ (0; 1). 1. Termination test. If z k satisÿes (1.1) – (1.2), then stop. 2. Direction. Find dk according to (1.3). 3. Step size. Let k be the largest step size for which z k + k dk ¿0. If k ¿1, set k = 1; otherwise set k =  k . 4. Update. Compute z k+1 = z k + k dk :
(1.4)
Increment k by one and return to Step 1.
2. Convergence To begin, we show that the sequence generated by the algorithm converges. Thereafter we prove, that the limit point solves (1.1) – (1.2) if a feasible solution exists. Several preliminary results are employed. The following three results characterize the update direction and convergence of the residual vector for (1.1). Lemma 1. The direction vector deÿned by (1:3) is given by dk = Z 2 QT (QZ 2 QT )−1 qk ;
(2.1)
k
k
k
Lemma 3. Let H ∈ Rn×(n−m) be a matrix whose columns form a basis for the null space of Q; and let h ∈ Rn such that Qh = q0 . Then dk = k [h − H (H T Z −2 H )−1 H T Z −2 h];
where dk is the direction given by (2:1); the scalars k are deÿned by (2:2); and Z = diag(z k ). Proof. By Lemma 2, qk =q−Qz k = k q0 = k Qh so that Q(z k + k h)=q. Hence, by deÿnition of H , any solution z for Qz = q is given by z = z k + k h − Hv so that dk = k h−Hv, for some v ∈ Rn−m . To determine v, problem (1.3) can be stated as minv∈Rn−m (h k −Hv)T Z −2 (h k − Hv). Optimization yields v= k (H T Z −2 H )−1 H T Z −2 h implying (2.3). Lemma 5 below implies that the set of all possible direction vectors is bounded. To prove this result we employ a preliminary result and the following notation. Let B = {{j1 ; : : : ; jm } | 16j1 ¡ · · · ¡ jm 6n}. let QB ∈ Rm×m be the submatrix deFor B ∈ B, ÿned by columns (indexed by) B of Q. Similarly, for I ∈ Rn×n , let IB ∈ Rn×m be the submatrix deÿned by columns B of I . Denote B = {B ∈ B | rank(QB ) = m}. Lemma 4. Let Q ∈ Rm×n ; rank(Q) = m and S = diag(s); s ∈ Rn ; s ¿ 0. Then X B IB QB−1 ; (2.4) SQT (QSQT )−1 = B∈B
where B = B =;
where q is the residual vector deÿned by q =q−Qz ; and Z = diag(z k ).
B = |QB |2
Proof. By straightforward optimization in (1.3).
and
Lemma 2. There exists a monotonically decreasing such that sequence { k } of scalars converging to ¿0 the residual vectors qk ; k = 1; 2; : : : ; satisfy k
k
q = q − Qz = (1 −
k−1
k−1
)q
k 0
= q :
(2.2)
Proof. Employing (1.4) – (2.1), qk+1 = q − Qz k+1 = qk − k Qdk = (1 − k )qk . Deÿne 0 = 1 and k+1 = (1 − k ) k , for k¿0. Observing that 0 ¡ k 61, the sequence { k } is nonnegative and monotonically decreasing and hence it converges to some ¿0.
(2.3)
(2.5) Y
sj
(2.6)
j∈B
= |QSQT | =
X B∈B
B =
X
B :
(2.7)
B∈B
replacing components Proof. For s ∈ Rn and B ∈ B, sj ; j 6∈ B, by zero, we obtain a vector denoted by sB . Let SB =diag(sB ) and SB =IBT SB IB (the diagonal matrix of elements sj ; j ∈ B). Then = (s), the determinant of QSQT , is a polynomial of s with the following properties: (i) each term of (s) is of degree m because each element of QSQT is a linear function of s; (sB ) = |QSB QT | = |QB SB QBT | = (ii) for all B ∈ B;
M. Kallio, S. Salo / Operations Research Letters 27 (2000) 101–107
|QB |2 |SB | = B ; (iii) if s has less than m nonzero components, then rank(QS) Q¡ m mand (s) = 0. Consider a term a j∈B sj j of (s) and suppose mj ¿ 1, for some j. Then (i) implies |B| ¡ m. Hence (iii) implies (s) = 0, for all s such that sj = 0, for j 6∈ B. It is elementary to show, that the polynomial (S) with such a Q property must have a = 0. Consequently, m in a term a j∈B sj j ; a 6= 0 implies mj = 1; j ∈ B, and P furthermore, B ∈ B by (i). Hence by (ii), (s) = B∈B B . But (sB ) = B = 0, for B ∈ B \ B, so that (2.7) follows. To prove (2.4), we employ the following known facts: if A and B are square matrices of equal size and A∗ and B∗ denote their adjoint matrices, then AA∗ = |A|I and (AB)∗ = B∗ A∗ . Denote P = SQT (QSQT )∗ . Elements of (QSQT )∗ are minors of order m − 1 of |QSQT |. Therefore, each element of the matrix P(s), is a polynomial of s and each term of such polynomials is of degree m. Using the facts above, for all P(sB )=SB QT (QSB QT )∗ =IB SB QBT (QB SB QBT )∗ = B ∈ B; |S B kQB |IB QB∗ . Hence, P(s)=0 for all s with s=sB ; B ∈ B \ B, or with s having less than m nonzero components. Repeating the same arguments as with (s), properties of P(s), with QB∗ = |QB |QB−1 , for B ∈ B, imply X B IB QB−1 ; (2.8) P= B∈B
where B is given by (2.6). Finally, observing that SQT (QSQT )−1 = −1 P, (2.8) yields (2.4). Corollary 1. For B ∈ B; let dkB be the basic solution to Qd = qk with dj ; j ∈ B as basic variables; i.e. dkB =IB QB−1 qk . Then the direction vector dk in Lemma 1 is the convex combination X B dkB ; (2.9) dk = B∈B
where B ’s are deÿned by (2:5)–(2:7) and S =Z 2 ; with Z = diag(z k ). Proof. By direct application of Lemma 4. Lemma 5. Let H and h be deÿned as in Lemma 3; let s ∈ Rn and deÿne S = diag(s) so that Se = s. For 0 ¡ s ¡ ∞; deÿne ds = h − H (H T S −2 H )−1 H T S −2 h. Then there is a scalar ¡ ∞ such that sups kds k ¡ .
103
2 T 2 T −1 0 Proof. By Lemmas 1 and 3, dP s = S Q (QS Q ) q 0 and thus by Corollary 1 ds = B∈B B dB , where the coecients B ∈ (0; 1)Pare deÿned by (2.5), B = Q |QB |2 j∈B sj2 and = B∈B B . Basic solutions d0B are bounded for all B ∈ B and the number of basic solutions is ÿnite. Thus kds k is bounded.
Theorem 1. The sequence {z k } generated by the algorithm converges. Proof. If the sequence {z k } is ÿnite, the algorithm stops in a solution z satisfying (1.1) – (1.2). Assume next that the sequence {z k } is inÿnite. By Lemma 5, there is 0 ¡ ∞ such that in (2.3), dk = k (h − r) for some r such that Qr = 0 and krk60 . Deÿne a closed cone of directions D = {d ∈ Rn |d = (h − r); ¿0; Qr = 0; krk60 } so that dk ∈ D for all k. By Lemma 2, the residuals qk converge to q0 . De For z ∈ Rn , ÿne M = {z ∈ Rn | q − Qz = q0 ; ¿ }. n denote z + D = {w ∈ R | w = z + d; d ∈ D}. Let C k = {z k + D} ∩ M so that z k ∈ C k . For some d ∈ D; z k+1 = z k + d ∈ C k , because z k+1 ∈ M . Because D is a cone, we have C k+1 = {z k+1 + D} ∩ M = {z k + d + D} ∩ M ⊂{z k + D} ∩ M = C k . Consequently, z l ∈ C l ⊂ C k , for l¿k. If z ∈ C k , then there is r ∈ Rn , with Qr = 0 and krk60 , and ¿0 such that z = z k + (h − r). Hence, employing Lemmas 2 and 3, q−Qz =q−Q(z k +(h− Thus, k − ¿¿0. r)) = ( k − )q0 , with k − ¿ . k By Lemma 2, − converges to 0, and therefore, the diameter of C k converges to 0. Hence, {z k } is a Cauchy sequence, wherefore it converges. Theorem 2. If a feasible solution zˆ for (1:1)–(1:2) exists; then {z k } converges to a feasible solution z. Furthermore; if {z k } is inÿnite; then zj = 0 implies ˆ zˆj = 0 for all feasible solutions z. Proof. Assume that a feasible solution zˆ exists. By If n = m, Lemma 1 Theorem 1, z k converges to z. ˆ If convergence is ÿnite, z is implies z 1 = z = z¿0. feasible and there is nothing to prove. Otherwise, after possible permutation of the variables, partition z T = (z1T ; z2T ), with zi ∈ Rni ; i = 1; 2, such that z1 = 0 and z2 ¿ 0. Suppose n1 = 0; i.e. z ¿ 0. If k ¡ 1 for all k, then (1.4) yields z k+1 = z k +  k dk , with  ¿ 0 and lim inf k k k dk k ¿ 0. This contradicts convergence of z k . Hence, convergence is ÿnite and z ¿ 0 is feasible.
104
M. Kallio, S. Salo / Operations Research Letters 27 (2000) 101–107
From now on we assume that n ¿ m; n1 ¿ 0 and {z } is inÿnite. For all k and j (j = 1; : : : ; n), deÿne jk = supl¿k |zjl − zj |, which converge to zero. Consider a collection of diagonal scaling matrices S =diag(s) ∈ Rn×n , where s ∈ V k = {s ∈ Rn | s¿0; |sj − zj |6jk for all j}. Then z k ∈ V k and Z = diag(z k ) belongs to this collection. Because jk is nonincreasing, for all j; z l ∈ V l ⊂ V k , for all l¿k. For all k, observing (2.3), deÿne a closed cone of directions k
Dk = {d ∈ Rn | d = S 2 u; ¿0; S 2 u = h − Hv; H T u = 0; s ∈ V k } so that dl ∈ Dl ⊂ Dk , for l¿k. Deÿne C k =z k +Dk so that C k is a closed cone with vertex at z k . Hence z l ∈ C l ⊂ C k , for all l¿k. Because C k is closed, z ∈ C k . We aim to show that, if ¿ 0 or = 0 and zˆ1 6= 0, then z 6∈ C k , for some k, wherefore = 0 and zˆ1 = 0. We proceed with the following decomposition: 0 S1 d1 ; S= ; d= 0 S2 d2 where di ∈ Rni and Si ∈ Rni ×ni , for i = 1; 2. Pick an index 6n1 so that z = 0. For all k, there is l¿k, such that zl = l . Because z ∈ C l ; z = z l + d, for some d ∈ Dl . Hence, −d1 =z1l − z1 =z1l =−S12 u1 ¿ 0 with ¿ 0 and the diagonal of S1 strictly positive. Consequently, for k large enough, the diagonal of S is strictly positive and d = (h − Hv), with = l − and H T u = H T S −2 (h − Hv) = 0. In the latter equation, premultiplication by vT yields vT H T S −2 (h − Hv) = 0:
(2.10)
By assumption, Qzˆ = q so that Q(zˆ − z l ) = l q0 . Let h=(z−z ˆ l )= l so that Qh=q0 . Denote ÿ== l =1− = l l so that ÿ ∈ (0; 1]. Then d = z − z = (h − Hv) = ÿ(zˆ − z l ) − Hv. Solving for Hv yields Hv = z˜ − z;
(2.11)
where z˜ ≡ ÿzˆ + (1 − ÿ)z l , and l
(h − Hv) = z − z :
(2.12)
With (2.11) – (2.12) we have −2 vT H T S −2 (h−Hv)= (z− ˜ z) T S −2 (z l − z). Decomposing the latter expression and observing (2.10) yields z˜T1 S1−2 z1l + (z˜2 − z2 )T S2−2 (z2l − z2 ) = 0:
(2.13)
Let k and l¿k increase without limit. Then, if n1 ¡ n; z2l converges to z2 ¿ 0 and S2 converges to
diag(z2 ) so that (z˜2 − z2 )T S2−2 (z2l − z2 ) converges to zero. For the ÿrst term in (2.13), observing that zl = l ¿s ¿ 0, we have z˜T1 S1−2 z l ¿ z˜ zl =s2 ¿z˜ =l :
(2.14)
Suppose ¿ 0. Then ÿ converges to zero and z˜ =l in (2.14) converges to one. Therefore, for some k large enough and l¿k such that zl =l , (2.13) does not have a solution, implying z 6∈ C l . This is a contradiction, wherefore = 0. For the second part of the assertion, suppose = 0 and zˆ ¿ 0. Then = l ; ÿ = 1 and z˜ = zˆ so that z˜ =l = zˆ =l , which diverges as l converges to zero. Therefore, for some k and l¿k, (2.13) does not have a solution, implying z 6∈ C l . This is a contradiction, wherefore zˆ = 0. 3. Infeasibility Detection of infeasibility may take place in alternative ways. First, if (1.1) – (1.2) does not have a feasible solution, then by Theorem 1 the sequence {z k } converges to a point z. By Lemma 2, the corresponding residual vector is given by q − Qz = q0 , where q0 is the residual vector at the initial point z 0 and ¿ 0 is the limit point of the nonnegative and monotone decreasing sequence { k } deÿned by Lemma 2. Thus, an indication for infeasibility is that k converges to a strictly positive level. However, in a ÿnite number of iterations, it may be dicult to judge when to conclude infeasibility and stop. Alternatively, we may modify the method of Section 1 (in a way to be discussed shortly) to ÿnd z ∈ Rn and Ä ∈ R to satisfy the homogeneous system Qz − Äq = 0;
(3.1)
z¿0; Ä¿0;
(3.2)
which is always feasible. Direct application of the method in Section 1 may not work as illustrated by the following example. Let (1.1) be given by z1 − 2z2 = 5 so that zˆ1 = 5; zˆ2 = 0 and Ĉ = 1 solve (3.1) – (3.2). Initializing our method for (3.1) – (3.2) with z10 = z20 = Ä0 = 1 yields in one iteration z11 =1:2; z21 =0:6 and Ä1 =0, which is feasible for (3.1) – (3.2) so that the iterations stop. However, this solution provides a homogeneous solution, but not a feasible solution for (1.1) – (1.2).
M. Kallio, S. Salo / Operations Research Letters 27 (2000) 101–107
To prevent such an early termination, we modify the step size rule of Step 3 in Section 1 by setting k = 1 if k ¿ 1 (instead of the earlier condition k ¿1). 1 Otherwise as before, k =  k , where  is the step size parameter. Then all results of Section 2 remain valid, but convergence is ÿnite only if a feasible interior point is encountered. Hence by Theorem 2, the modiÿed algorithm applied to (3.1) – (3.2) converges to a feasible solution z and Ä. If the original system (1.1) – (1.2) has a feasible solution z, ˆ then zˆ and Ĉ = 1 solve (3.1) – (3.2). In such a case, if convergence is ÿnite, we obtain an interior solution z ¿ 0 and Ä ¿ 0, and if convergence is inÿnite, Ĉ ¿ 0 implies Ä ¿ 0 by Theorem 2. Consequently in both cases, z= Ä solves (1.1) – (1.2). In the example above, with  = 0:99, ÿnite convergence occurs in iteration two, and rounded ÿgures are z1 = 1:228; z2 = 0:589 and Ä = 0:00999. Conversely, if (1.1) – (1.2) is infeasible, then the solution z and Ä for the homogeneous system (3.1) – (3.2) obtained by the method of Section 1, employing the modiÿed step size rule, satisÿes Ä = 0. 4. Free and bounded variables We consider next a modiÿcation of the system (1.1) – (1.2), where some of the variables are free and some have simple bounds. A free variable zj has an upper bound equal to ∞ and a lower bound equal to −∞. If a variable zj is nonfree, but its lower bound is −∞, we replace zj by −zj , so that the revised variable has a ÿnite lower bound. Thereafter, we translate each nonfree variable so that its lower bound becomes 0; i.e. if a variable zj has a lower bound lj ¿ − ∞, we replace zj by zj + lj . Thereby, if the revised upper bound is less than ∞, it is denoted by uj . Finally, we subdivide the set of indices J = {1; 2; : : : ; n} into J = JF ∪ JB ∪ JP so that −∞6zj 6∞ 06zj 6uj 06zj 6∞
for j ∈ JF ; for j ∈ JB ; for j ∈ JP :
Hence, in the resulting modiÿed system, (1.2) is replaced by 06zj 6uj , for j ∈ JB , and 06zj , for j ∈ JP . 1 In an alternative modiÿcation we set k = 1 if k ¿1, provided that the component Ä does not decrease to zero in iteration k. To shorten the discussion, we do not pursue this alternative further.
105
Introducing new variables in the usual way, an equivalent system can be reformulated employing nonnegative variables only. For free variables j ∈ JF , substitute zj by zj+ − zj− with zj+ ¿0 and zj− ¿0. For bounded variables j ∈ JB , introduce a slack variable vj so that the requirement zj 6uj becomes zj + vj = uj with vj ¿0. Denoting by Qj the column vector j of Q, the system is now revised to X X X Qj (zj+ − zj− ) + Qj zj + Qj zj = q; (4.1) j∈JF
j∈JB
zj + vj = uj
for j ∈ JB ;
zj+ ; zj− ¿0
for j ∈ JF ;
j∈JP
(4.2) zj ; vj ¿0 for j ∈ JB ;
zj ¿0 for j ∈ Jp :
(4.3)
Clearly, the method of Section 1 applies to solve system (4.1) – (4.3). Straightforward algebra shows, that the following steps can be employed to compute the direction vector of update in each iteration k. Suppressing the iteration count k for variables zj+ ; zj− ; vj and zj , deÿne a diagonal scaling matrix S = diag(s) ∈ Rn×n as follows: + 2 − 2 for j ∈ JF ; (zj ) + (zj ) 2 2 2 sj = (vj zj ) =(vj + zj ) for j ∈ JB ; 2 for j ∈ JP : zj Deÿne residuals, ukj = uj − zj − vj , for j ∈ JB ; and P qˆk = qk − j∈JB (sj =vj2 )Qj ukj . Deÿne an auxiliary vector k ∈ Rm by k = (QSQT )−1 qˆk and scalars jk , for j ∈ JB , by jk = (u j − zj2 QjT k )=(zj2 + vj2 ). Then the components of the direction vector, obtained by applying (2.1) to the revised system (4.1) – (4.3), are as follows: (zj+ )2 QjT k for zj+ ; j ∈ JF ; − 2 T k for zj− ; j ∈ JF ; −(zj ) Qj dkj =
zj2 (QjT k + jk ) vj2 jk z 2 Q T k j j
for zj ; j ∈ JB ;
for vj ; j ∈ JB ; for zj ; j ∈ JP :
Note that the essential part of computational e ort relates to the inverse of the matrix QSQT , which has the same dimension and structure as the matrix QZ 2 QT in (2.1) of the original system (1.1) – (1.2).
106
M. Kallio, S. Salo / Operations Research Letters 27 (2000) 101–107 Table 1 The number of rows (Rows), columns (Cols) of test problems, as well as the number of iterations for Andersen (AA) and the method (KS) of Section 1 Iterations Problem
Rows
Cols
AA
KS
aÿro sc50b sc50a sc105 adlittle scagr7 stocfor1 blend sc205 share2b lotÿ share1b scorpion scagr25 sctap1 brandy israel scsd1 agg bandm e226 scfxm1 beaconfd degen2 agg2 agg3 scrs8 scsd6 ship04s bnl1 scfxm2
27 50 50 105 56 129 117 74 295 96 153 117 388 471 300 220 174 77 488 305 223 330 173 444 516 516 490 147 402 643 660
32 48 48 103 97 140 111 83 203 79 308 225 358 500 480 249 142 760 163 472 282 457 262 534 302 302 1169 1350 1458 1175 914
9 11 11 12 13 15 12 10 12 11 17 26 14 16 19 14 19 10 17 15 16 18 12 10 18 17 25 12 17 26 20
21 18 20 21 24 27 24 24 23 27 30 38 21 40 32 30 48 25 72 30 35 38 25 27 34 35 59 27 25 78 41
5. Application to linear programming Let A ∈ Rk×l ; c ∈ Rl and b ∈ Rk , and consider the following linear programming problem: min{cT x | Ax = b; x¿0} x∈Rl
(5.1)
and its dual max {bT y | AT y + s = c; s¿0}:
y∈Rk ; s∈R1
(5.2)
Optimality conditions for (5.1) – (5.2) may be stated by (1.1) with
0 Q= A −cT
AT 0 bT
I 0; 0
c q = b; 0
x z=y s
and by nonnegativity constraints x¿0 and s¿0. To deal with the free variables y and to solve the problem, we may proceed as described in Section 4. A rudimentary code is set up to test the method for some linear programming test problems from the Netlib collection. With this code, we only aim to study the number of iterations needed to solve the problem. Otherwise we do not worry about computational eciency and memory requirements at this stage. Consequently, we employ dense Cholesky factorization to deal with the inverse (QSQT )−1 and we use high
M. Kallio, S. Salo / Operations Research Letters 27 (2000) 101–107
enough accuracy to deal with ill conditioning of the matrix QSQT . We prescale the matrix A with a procedure that aims to equalize the norms of all rows and columns of A close to 1. Thereafter, to initialize the iterations, we set all variables equal to ten. Our experimental code does not deal with simple bounds; i.e. all variables are nonnegative or free. We compare the results with a state of the art interior point code by Andersen and Andersen [1], wherefore a similar stopping criterion is employed in both codes. In particular, the vital requirement is that the absolute value of the relative duality gap |cT x −bT y|=(1+|cT x|) is less than 10−8 . Besides, both codes employ primal and dual infeasibility tolerances, which di er, because initialization procedures di er. We require the norms of primal and dual infeasibility vectors relative to norms of b and c; respectively, to be less than 10−6 . For the step size parameter, we use  = 0:99. The results are shown in Table 1 below. We obtain approximately twice as many iterations as compared with using MOSEK v1.0.beta [1]. Given that the latter code employs most known enhancements for primal–
107
dual infeasible interior point methods, we may regard our results quite satisfactory. Acknowledgements The authors wish to thank Erling Andersen for providing them with the results obtained with MOSEK, as well as an anonymous referee for drawing their attention to infeasible problems. Financial support for this work is gratefully acknowledged from The Foundation for the Helsinki School of Economics. References [1] E.D. Andersen, K.D. Andersen, The MOSEK interior point optimizer for linear programming: an implementation of the homogeneous algorithm, in: H. Frenk, K. Roos, T. Terlaky, S. Zhang (Eds.), High Performance Optimization, Proceedings of the HPOPT-II Conference, 2000. [2] T. Terlaky (Ed.), Interior Point Methods of Mathematical Programming, Kluwer Academic Publishers, Dordrecht, 1996.
www.elsevier.com/locate/dsw
An interior point method for solving systems of linear equations and inequalities Markku Kallio ∗ , Seppo Salo Helsinki School of Economics, Runeberginkatu 14-16, Helsinki 00100, Finland Received 1 March 1999; received in revised form 1 March 2000
Abstract A simple interior point method is proposed for solving a system of linear equations subject to nonnegativity constraints. The direction of update is deÿned by projection of the current solution on a linear manifold deÿned by the equations. Infeasibility is discussed and extension for free and bounded variables is presented. As an application, we consider linear c programming problems and a comparison with a state-of-the-art primal–dual infeasible interior point code is presented. 2000 Elsevier Science B.V. All rights reserved. Keywords: Interior point algorithms; Linear inequalities; Linear programming
1. The problem and a method of solution An inspiration to our article arises from the vast literature on interior point methods for solving mathematical programming problems (see e.g. [2]). In particular, we are motivated by the success of primal–dual infeasible interior point algorithms for linear programming. Our aim is to solve systems of linear equations subject to nonnegativity constraints, for which optimality conditions for linear programming is a special case. We propose a simple interior point algorithm which can be initiated with any positive vector. It always converges, and if the problem is feasible, the limit point is a solution.
Let q ∈ Rm and Q ∈ Rm×n with m ¡ n and rank(Q) = m. Our problem is to ÿnd z ∈ Rn such that Qz = q; z¿0:
Corresponding author. E-mail address: kallio@hkkk.ÿ (M. Kallio).
(1.2) n
Denote M = {z ∈ R | Qz = q}; Z = diag(z), and for d ∈ Rn , deÿne Z −2 -norm by kdkZ −2 = (dT Z −2 d)1=2 . Our method starts with z 0 ∈ Rn ; z 0 ¿ 0. For iteration k; k = 0; 1; 2; : : : ; denote Z = diag(z k ), suppressing k in Z. The direction dk of update is deÿned by dk = z − z k where z is the point of M which is closest to z k in Z −2 -norm; i.e., z is the unique solution to the following problem: min kz − z k kZ −2 : z∈M
∗
(1.1)
(1.3)
A step in this direction is taken. If z k + dk ¿0 a solution for (1.1) – (1.2) is found; otherwise the step size k ¡ 1 is applied so that z k+1 = z k + k dk ¿ 0.
c 2000 Elsevier Science B.V. All rights reserved. 0167-6377/00/$ - see front matter PII: S 0 1 6 7 - 6 3 7 7 ( 0 0 ) 0 0 0 3 6 - 5
102
M. Kallio, S. Salo / Operations Research Letters 27 (2000) 101–107
The steps of the algorithm are formalized as follows: 0. Initialization. Set the iteration count k to 0, pick any z 0 ∈ Rn such that z 0 ¿ 0, and let the step size parameter  ∈ (0; 1). 1. Termination test. If z k satisÿes (1.1) – (1.2), then stop. 2. Direction. Find dk according to (1.3). 3. Step size. Let k be the largest step size for which z k + k dk ¿0. If k ¿1, set k = 1; otherwise set k =  k . 4. Update. Compute z k+1 = z k + k dk :
(1.4)
Increment k by one and return to Step 1.
2. Convergence To begin, we show that the sequence generated by the algorithm converges. Thereafter we prove, that the limit point solves (1.1) – (1.2) if a feasible solution exists. Several preliminary results are employed. The following three results characterize the update direction and convergence of the residual vector for (1.1). Lemma 1. The direction vector deÿned by (1:3) is given by dk = Z 2 QT (QZ 2 QT )−1 qk ;
(2.1)
k
k
k
Lemma 3. Let H ∈ Rn×(n−m) be a matrix whose columns form a basis for the null space of Q; and let h ∈ Rn such that Qh = q0 . Then dk = k [h − H (H T Z −2 H )−1 H T Z −2 h];
where dk is the direction given by (2:1); the scalars k are deÿned by (2:2); and Z = diag(z k ). Proof. By Lemma 2, qk =q−Qz k = k q0 = k Qh so that Q(z k + k h)=q. Hence, by deÿnition of H , any solution z for Qz = q is given by z = z k + k h − Hv so that dk = k h−Hv, for some v ∈ Rn−m . To determine v, problem (1.3) can be stated as minv∈Rn−m (h k −Hv)T Z −2 (h k − Hv). Optimization yields v= k (H T Z −2 H )−1 H T Z −2 h implying (2.3). Lemma 5 below implies that the set of all possible direction vectors is bounded. To prove this result we employ a preliminary result and the following notation. Let B = {{j1 ; : : : ; jm } | 16j1 ¡ · · · ¡ jm 6n}. let QB ∈ Rm×m be the submatrix deFor B ∈ B, ÿned by columns (indexed by) B of Q. Similarly, for I ∈ Rn×n , let IB ∈ Rn×m be the submatrix deÿned by columns B of I . Denote B = {B ∈ B | rank(QB ) = m}. Lemma 4. Let Q ∈ Rm×n ; rank(Q) = m and S = diag(s); s ∈ Rn ; s ¿ 0. Then X B IB QB−1 ; (2.4) SQT (QSQT )−1 = B∈B
where B = B =;
where q is the residual vector deÿned by q =q−Qz ; and Z = diag(z k ).
B = |QB |2
Proof. By straightforward optimization in (1.3).
and
Lemma 2. There exists a monotonically decreasing such that sequence { k } of scalars converging to ¿0 the residual vectors qk ; k = 1; 2; : : : ; satisfy k
k
q = q − Qz = (1 −
k−1
k−1
)q
k 0
= q :
(2.2)
Proof. Employing (1.4) – (2.1), qk+1 = q − Qz k+1 = qk − k Qdk = (1 − k )qk . Deÿne 0 = 1 and k+1 = (1 − k ) k , for k¿0. Observing that 0 ¡ k 61, the sequence { k } is nonnegative and monotonically decreasing and hence it converges to some ¿0.
(2.3)
(2.5) Y
sj
(2.6)
j∈B
= |QSQT | =
X B∈B
B =
X
B :
(2.7)
B∈B
replacing components Proof. For s ∈ Rn and B ∈ B, sj ; j 6∈ B, by zero, we obtain a vector denoted by sB . Let SB =diag(sB ) and SB =IBT SB IB (the diagonal matrix of elements sj ; j ∈ B). Then = (s), the determinant of QSQT , is a polynomial of s with the following properties: (i) each term of (s) is of degree m because each element of QSQT is a linear function of s; (sB ) = |QSB QT | = |QB SB QBT | = (ii) for all B ∈ B;
M. Kallio, S. Salo / Operations Research Letters 27 (2000) 101–107
|QB |2 |SB | = B ; (iii) if s has less than m nonzero components, then rank(QS) Q¡ m mand (s) = 0. Consider a term a j∈B sj j of (s) and suppose mj ¿ 1, for some j. Then (i) implies |B| ¡ m. Hence (iii) implies (s) = 0, for all s such that sj = 0, for j 6∈ B. It is elementary to show, that the polynomial (S) with such a Q property must have a = 0. Consequently, m in a term a j∈B sj j ; a 6= 0 implies mj = 1; j ∈ B, and P furthermore, B ∈ B by (i). Hence by (ii), (s) = B∈B B . But (sB ) = B = 0, for B ∈ B \ B, so that (2.7) follows. To prove (2.4), we employ the following known facts: if A and B are square matrices of equal size and A∗ and B∗ denote their adjoint matrices, then AA∗ = |A|I and (AB)∗ = B∗ A∗ . Denote P = SQT (QSQT )∗ . Elements of (QSQT )∗ are minors of order m − 1 of |QSQT |. Therefore, each element of the matrix P(s), is a polynomial of s and each term of such polynomials is of degree m. Using the facts above, for all P(sB )=SB QT (QSB QT )∗ =IB SB QBT (QB SB QBT )∗ = B ∈ B; |S B kQB |IB QB∗ . Hence, P(s)=0 for all s with s=sB ; B ∈ B \ B, or with s having less than m nonzero components. Repeating the same arguments as with (s), properties of P(s), with QB∗ = |QB |QB−1 , for B ∈ B, imply X B IB QB−1 ; (2.8) P= B∈B
where B is given by (2.6). Finally, observing that SQT (QSQT )−1 = −1 P, (2.8) yields (2.4). Corollary 1. For B ∈ B; let dkB be the basic solution to Qd = qk with dj ; j ∈ B as basic variables; i.e. dkB =IB QB−1 qk . Then the direction vector dk in Lemma 1 is the convex combination X B dkB ; (2.9) dk = B∈B
where B ’s are deÿned by (2:5)–(2:7) and S =Z 2 ; with Z = diag(z k ). Proof. By direct application of Lemma 4. Lemma 5. Let H and h be deÿned as in Lemma 3; let s ∈ Rn and deÿne S = diag(s) so that Se = s. For 0 ¡ s ¡ ∞; deÿne ds = h − H (H T S −2 H )−1 H T S −2 h. Then there is a scalar ¡ ∞ such that sups kds k ¡ .
103
2 T 2 T −1 0 Proof. By Lemmas 1 and 3, dP s = S Q (QS Q ) q 0 and thus by Corollary 1 ds = B∈B B dB , where the coecients B ∈ (0; 1)Pare deÿned by (2.5), B = Q |QB |2 j∈B sj2 and = B∈B B . Basic solutions d0B are bounded for all B ∈ B and the number of basic solutions is ÿnite. Thus kds k is bounded.
Theorem 1. The sequence {z k } generated by the algorithm converges. Proof. If the sequence {z k } is ÿnite, the algorithm stops in a solution z satisfying (1.1) – (1.2). Assume next that the sequence {z k } is inÿnite. By Lemma 5, there is 0 ¡ ∞ such that in (2.3), dk = k (h − r) for some r such that Qr = 0 and krk60 . Deÿne a closed cone of directions D = {d ∈ Rn |d = (h − r); ¿0; Qr = 0; krk60 } so that dk ∈ D for all k. By Lemma 2, the residuals qk converge to q0 . De For z ∈ Rn , ÿne M = {z ∈ Rn | q − Qz = q0 ; ¿ }. n denote z + D = {w ∈ R | w = z + d; d ∈ D}. Let C k = {z k + D} ∩ M so that z k ∈ C k . For some d ∈ D; z k+1 = z k + d ∈ C k , because z k+1 ∈ M . Because D is a cone, we have C k+1 = {z k+1 + D} ∩ M = {z k + d + D} ∩ M ⊂{z k + D} ∩ M = C k . Consequently, z l ∈ C l ⊂ C k , for l¿k. If z ∈ C k , then there is r ∈ Rn , with Qr = 0 and krk60 , and ¿0 such that z = z k + (h − r). Hence, employing Lemmas 2 and 3, q−Qz =q−Q(z k +(h− Thus, k − ¿¿0. r)) = ( k − )q0 , with k − ¿ . k By Lemma 2, − converges to 0, and therefore, the diameter of C k converges to 0. Hence, {z k } is a Cauchy sequence, wherefore it converges. Theorem 2. If a feasible solution zˆ for (1:1)–(1:2) exists; then {z k } converges to a feasible solution z. Furthermore; if {z k } is inÿnite; then zj = 0 implies ˆ zˆj = 0 for all feasible solutions z. Proof. Assume that a feasible solution zˆ exists. By If n = m, Lemma 1 Theorem 1, z k converges to z. ˆ If convergence is ÿnite, z is implies z 1 = z = z¿0. feasible and there is nothing to prove. Otherwise, after possible permutation of the variables, partition z T = (z1T ; z2T ), with zi ∈ Rni ; i = 1; 2, such that z1 = 0 and z2 ¿ 0. Suppose n1 = 0; i.e. z ¿ 0. If k ¡ 1 for all k, then (1.4) yields z k+1 = z k +  k dk , with  ¿ 0 and lim inf k k k dk k ¿ 0. This contradicts convergence of z k . Hence, convergence is ÿnite and z ¿ 0 is feasible.
104
M. Kallio, S. Salo / Operations Research Letters 27 (2000) 101–107
From now on we assume that n ¿ m; n1 ¿ 0 and {z } is inÿnite. For all k and j (j = 1; : : : ; n), deÿne jk = supl¿k |zjl − zj |, which converge to zero. Consider a collection of diagonal scaling matrices S =diag(s) ∈ Rn×n , where s ∈ V k = {s ∈ Rn | s¿0; |sj − zj |6jk for all j}. Then z k ∈ V k and Z = diag(z k ) belongs to this collection. Because jk is nonincreasing, for all j; z l ∈ V l ⊂ V k , for all l¿k. For all k, observing (2.3), deÿne a closed cone of directions k
Dk = {d ∈ Rn | d = S 2 u; ¿0; S 2 u = h − Hv; H T u = 0; s ∈ V k } so that dl ∈ Dl ⊂ Dk , for l¿k. Deÿne C k =z k +Dk so that C k is a closed cone with vertex at z k . Hence z l ∈ C l ⊂ C k , for all l¿k. Because C k is closed, z ∈ C k . We aim to show that, if ¿ 0 or = 0 and zˆ1 6= 0, then z 6∈ C k , for some k, wherefore = 0 and zˆ1 = 0. We proceed with the following decomposition: 0 S1 d1 ; S= ; d= 0 S2 d2 where di ∈ Rni and Si ∈ Rni ×ni , for i = 1; 2. Pick an index 6n1 so that z = 0. For all k, there is l¿k, such that zl = l . Because z ∈ C l ; z = z l + d, for some d ∈ Dl . Hence, −d1 =z1l − z1 =z1l =−S12 u1 ¿ 0 with ¿ 0 and the diagonal of S1 strictly positive. Consequently, for k large enough, the diagonal of S is strictly positive and d = (h − Hv), with = l − and H T u = H T S −2 (h − Hv) = 0. In the latter equation, premultiplication by vT yields vT H T S −2 (h − Hv) = 0:
(2.10)
By assumption, Qzˆ = q so that Q(zˆ − z l ) = l q0 . Let h=(z−z ˆ l )= l so that Qh=q0 . Denote ÿ== l =1− = l l so that ÿ ∈ (0; 1]. Then d = z − z = (h − Hv) = ÿ(zˆ − z l ) − Hv. Solving for Hv yields Hv = z˜ − z;
(2.11)
where z˜ ≡ ÿzˆ + (1 − ÿ)z l , and l
(h − Hv) = z − z :
(2.12)
With (2.11) – (2.12) we have −2 vT H T S −2 (h−Hv)= (z− ˜ z) T S −2 (z l − z). Decomposing the latter expression and observing (2.10) yields z˜T1 S1−2 z1l + (z˜2 − z2 )T S2−2 (z2l − z2 ) = 0:
(2.13)
Let k and l¿k increase without limit. Then, if n1 ¡ n; z2l converges to z2 ¿ 0 and S2 converges to
diag(z2 ) so that (z˜2 − z2 )T S2−2 (z2l − z2 ) converges to zero. For the ÿrst term in (2.13), observing that zl = l ¿s ¿ 0, we have z˜T1 S1−2 z l ¿ z˜ zl =s2 ¿z˜ =l :
(2.14)
Suppose ¿ 0. Then ÿ converges to zero and z˜ =l in (2.14) converges to one. Therefore, for some k large enough and l¿k such that zl =l , (2.13) does not have a solution, implying z 6∈ C l . This is a contradiction, wherefore = 0. For the second part of the assertion, suppose = 0 and zˆ ¿ 0. Then = l ; ÿ = 1 and z˜ = zˆ so that z˜ =l = zˆ =l , which diverges as l converges to zero. Therefore, for some k and l¿k, (2.13) does not have a solution, implying z 6∈ C l . This is a contradiction, wherefore zˆ = 0. 3. Infeasibility Detection of infeasibility may take place in alternative ways. First, if (1.1) – (1.2) does not have a feasible solution, then by Theorem 1 the sequence {z k } converges to a point z. By Lemma 2, the corresponding residual vector is given by q − Qz = q0 , where q0 is the residual vector at the initial point z 0 and ¿ 0 is the limit point of the nonnegative and monotone decreasing sequence { k } deÿned by Lemma 2. Thus, an indication for infeasibility is that k converges to a strictly positive level. However, in a ÿnite number of iterations, it may be dicult to judge when to conclude infeasibility and stop. Alternatively, we may modify the method of Section 1 (in a way to be discussed shortly) to ÿnd z ∈ Rn and Ä ∈ R to satisfy the homogeneous system Qz − Äq = 0;
(3.1)
z¿0; Ä¿0;
(3.2)
which is always feasible. Direct application of the method in Section 1 may not work as illustrated by the following example. Let (1.1) be given by z1 − 2z2 = 5 so that zˆ1 = 5; zˆ2 = 0 and Ĉ = 1 solve (3.1) – (3.2). Initializing our method for (3.1) – (3.2) with z10 = z20 = Ä0 = 1 yields in one iteration z11 =1:2; z21 =0:6 and Ä1 =0, which is feasible for (3.1) – (3.2) so that the iterations stop. However, this solution provides a homogeneous solution, but not a feasible solution for (1.1) – (1.2).
M. Kallio, S. Salo / Operations Research Letters 27 (2000) 101–107
To prevent such an early termination, we modify the step size rule of Step 3 in Section 1 by setting k = 1 if k ¿ 1 (instead of the earlier condition k ¿1). 1 Otherwise as before, k =  k , where  is the step size parameter. Then all results of Section 2 remain valid, but convergence is ÿnite only if a feasible interior point is encountered. Hence by Theorem 2, the modiÿed algorithm applied to (3.1) – (3.2) converges to a feasible solution z and Ä. If the original system (1.1) – (1.2) has a feasible solution z, ˆ then zˆ and Ĉ = 1 solve (3.1) – (3.2). In such a case, if convergence is ÿnite, we obtain an interior solution z ¿ 0 and Ä ¿ 0, and if convergence is inÿnite, Ĉ ¿ 0 implies Ä ¿ 0 by Theorem 2. Consequently in both cases, z= Ä solves (1.1) – (1.2). In the example above, with  = 0:99, ÿnite convergence occurs in iteration two, and rounded ÿgures are z1 = 1:228; z2 = 0:589 and Ä = 0:00999. Conversely, if (1.1) – (1.2) is infeasible, then the solution z and Ä for the homogeneous system (3.1) – (3.2) obtained by the method of Section 1, employing the modiÿed step size rule, satisÿes Ä = 0. 4. Free and bounded variables We consider next a modiÿcation of the system (1.1) – (1.2), where some of the variables are free and some have simple bounds. A free variable zj has an upper bound equal to ∞ and a lower bound equal to −∞. If a variable zj is nonfree, but its lower bound is −∞, we replace zj by −zj , so that the revised variable has a ÿnite lower bound. Thereafter, we translate each nonfree variable so that its lower bound becomes 0; i.e. if a variable zj has a lower bound lj ¿ − ∞, we replace zj by zj + lj . Thereby, if the revised upper bound is less than ∞, it is denoted by uj . Finally, we subdivide the set of indices J = {1; 2; : : : ; n} into J = JF ∪ JB ∪ JP so that −∞6zj 6∞ 06zj 6uj 06zj 6∞
for j ∈ JF ; for j ∈ JB ; for j ∈ JP :
Hence, in the resulting modiÿed system, (1.2) is replaced by 06zj 6uj , for j ∈ JB , and 06zj , for j ∈ JP . 1 In an alternative modiÿcation we set k = 1 if k ¿1, provided that the component Ä does not decrease to zero in iteration k. To shorten the discussion, we do not pursue this alternative further.
105
Introducing new variables in the usual way, an equivalent system can be reformulated employing nonnegative variables only. For free variables j ∈ JF , substitute zj by zj+ − zj− with zj+ ¿0 and zj− ¿0. For bounded variables j ∈ JB , introduce a slack variable vj so that the requirement zj 6uj becomes zj + vj = uj with vj ¿0. Denoting by Qj the column vector j of Q, the system is now revised to X X X Qj (zj+ − zj− ) + Qj zj + Qj zj = q; (4.1) j∈JF
j∈JB
zj + vj = uj
for j ∈ JB ;
zj+ ; zj− ¿0
for j ∈ JF ;
j∈JP
(4.2) zj ; vj ¿0 for j ∈ JB ;
zj ¿0 for j ∈ Jp :
(4.3)
Clearly, the method of Section 1 applies to solve system (4.1) – (4.3). Straightforward algebra shows, that the following steps can be employed to compute the direction vector of update in each iteration k. Suppressing the iteration count k for variables zj+ ; zj− ; vj and zj , deÿne a diagonal scaling matrix S = diag(s) ∈ Rn×n as follows: + 2 − 2 for j ∈ JF ; (zj ) + (zj ) 2 2 2 sj = (vj zj ) =(vj + zj ) for j ∈ JB ; 2 for j ∈ JP : zj Deÿne residuals, ukj = uj − zj − vj , for j ∈ JB ; and P qˆk = qk − j∈JB (sj =vj2 )Qj ukj . Deÿne an auxiliary vector k ∈ Rm by k = (QSQT )−1 qˆk and scalars jk , for j ∈ JB , by jk = (u j − zj2 QjT k )=(zj2 + vj2 ). Then the components of the direction vector, obtained by applying (2.1) to the revised system (4.1) – (4.3), are as follows: (zj+ )2 QjT k for zj+ ; j ∈ JF ; − 2 T k for zj− ; j ∈ JF ; −(zj ) Qj dkj =
zj2 (QjT k + jk ) vj2 jk z 2 Q T k j j
for zj ; j ∈ JB ;
for vj ; j ∈ JB ; for zj ; j ∈ JP :
Note that the essential part of computational e ort relates to the inverse of the matrix QSQT , which has the same dimension and structure as the matrix QZ 2 QT in (2.1) of the original system (1.1) – (1.2).
106
M. Kallio, S. Salo / Operations Research Letters 27 (2000) 101–107 Table 1 The number of rows (Rows), columns (Cols) of test problems, as well as the number of iterations for Andersen (AA) and the method (KS) of Section 1 Iterations Problem
Rows
Cols
AA
KS
aÿro sc50b sc50a sc105 adlittle scagr7 stocfor1 blend sc205 share2b lotÿ share1b scorpion scagr25 sctap1 brandy israel scsd1 agg bandm e226 scfxm1 beaconfd degen2 agg2 agg3 scrs8 scsd6 ship04s bnl1 scfxm2
27 50 50 105 56 129 117 74 295 96 153 117 388 471 300 220 174 77 488 305 223 330 173 444 516 516 490 147 402 643 660
32 48 48 103 97 140 111 83 203 79 308 225 358 500 480 249 142 760 163 472 282 457 262 534 302 302 1169 1350 1458 1175 914
9 11 11 12 13 15 12 10 12 11 17 26 14 16 19 14 19 10 17 15 16 18 12 10 18 17 25 12 17 26 20
21 18 20 21 24 27 24 24 23 27 30 38 21 40 32 30 48 25 72 30 35 38 25 27 34 35 59 27 25 78 41
5. Application to linear programming Let A ∈ Rk×l ; c ∈ Rl and b ∈ Rk , and consider the following linear programming problem: min{cT x | Ax = b; x¿0} x∈Rl
(5.1)
and its dual max {bT y | AT y + s = c; s¿0}:
y∈Rk ; s∈R1
(5.2)
Optimality conditions for (5.1) – (5.2) may be stated by (1.1) with
0 Q= A −cT
AT 0 bT
I 0; 0
c q = b; 0
x z=y s
and by nonnegativity constraints x¿0 and s¿0. To deal with the free variables y and to solve the problem, we may proceed as described in Section 4. A rudimentary code is set up to test the method for some linear programming test problems from the Netlib collection. With this code, we only aim to study the number of iterations needed to solve the problem. Otherwise we do not worry about computational eciency and memory requirements at this stage. Consequently, we employ dense Cholesky factorization to deal with the inverse (QSQT )−1 and we use high
M. Kallio, S. Salo / Operations Research Letters 27 (2000) 101–107
enough accuracy to deal with ill conditioning of the matrix QSQT . We prescale the matrix A with a procedure that aims to equalize the norms of all rows and columns of A close to 1. Thereafter, to initialize the iterations, we set all variables equal to ten. Our experimental code does not deal with simple bounds; i.e. all variables are nonnegative or free. We compare the results with a state of the art interior point code by Andersen and Andersen [1], wherefore a similar stopping criterion is employed in both codes. In particular, the vital requirement is that the absolute value of the relative duality gap |cT x −bT y|=(1+|cT x|) is less than 10−8 . Besides, both codes employ primal and dual infeasibility tolerances, which di er, because initialization procedures di er. We require the norms of primal and dual infeasibility vectors relative to norms of b and c; respectively, to be less than 10−6 . For the step size parameter, we use  = 0:99. The results are shown in Table 1 below. We obtain approximately twice as many iterations as compared with using MOSEK v1.0.beta [1]. Given that the latter code employs most known enhancements for primal–
107
dual infeasible interior point methods, we may regard our results quite satisfactory. Acknowledgements The authors wish to thank Erling Andersen for providing them with the results obtained with MOSEK, as well as an anonymous referee for drawing their attention to infeasible problems. Financial support for this work is gratefully acknowledged from The Foundation for the Helsinki School of Economics. References [1] E.D. Andersen, K.D. Andersen, The MOSEK interior point optimizer for linear programming: an implementation of the homogeneous algorithm, in: H. Frenk, K. Roos, T. Terlaky, S. Zhang (Eds.), High Performance Optimization, Proceedings of the HPOPT-II Conference, 2000. [2] T. Terlaky (Ed.), Interior Point Methods of Mathematical Programming, Kluwer Academic Publishers, Dordrecht, 1996.