Applied Mathematics and Computation 183 (2006) 646–658 www.elsevier.com/locate/amc
Polynomial time second order Mehrotra-type predictor–corrector algorithms Maziar Salahi a
a,*
, Nezam Mahdavi-Amiri
b
Department of Mathematics, Faculty of Sciences, The University of Guilan, Namjoo Street, Rasht, Iran b Department of Mathematical Sciences, Sharif University of Technology, Tehran, Iran
Abstract Salahi et al. [M. Salahi, J. Peng, T. Terlaky, On Mehrtora type predictor–corrector algorithms, Technical Report 2005/ 4, Advanced Optimization Lab, Department of Computing and Software, McMaster University, http://www.cas.mcmaster.ca/~oplab/publication, SIAM Journal on Optimization, submitted for publication] give a numerical example showing that Mehrotra’s original predictor–corrector algorithm, which is the basis of interior point methods software packages, may be very inefficient in practice. This motivated Salahi et al. to come up with a safeguarded algorithm that enjoys a polynomial iteration complexity and is efficient in practice. Here we discuss a variation of Mehrotra’s second order predictor– corrector algorithm [S. Mehrotra, On the implementation of a (primal–dual) interior point method, SIAM Journal on Optimization 2 (1992) 575–601] and use the example of Salahi et al. to show that the algorithm may have to take very small steps in order to remain in a certain neighborhood of the central path and subsequently needs excessively large number of iterations to stop. We then introduce a safeguard that guarantees a lower bound for the maximum step size in the corrector step of the algorithm and subsequently a polynomial number of iterations. A second modification of algorithm is proposed which enjoys even a better iteration complexity. Some limited encouraging computational results are reported. 2006 Elsevier Inc. All rights reserved. Keywords: Linear optimization; Predictor–corrector methods; Interior point methods; Second order methods; Polynomial complexity
1. Introduction Mehrotra-type predictor–corrector algorithms are the base of the interior point methods (IPMs) software packages such as [1,11,12] and many others. Due to the practical efficiency of the algorithm the authors of [7] analyzed a feasible variant of Mehrotra’s original algorithm. By a numerical example they showed that this algorithm which is using an adaptive update of the barrier (centering) parameter at each iteration might be very inefficient in practice. This observation motivated them to combine the algorithm with some safeguards to prevent such a phenomenon. In this paper we analyze a feasible version of a variation of Mehrotra’s second * Corresponding author. Address: Advanced Optimization Laboratory, McMaster University, 1280 Main Street West, Hamilton, Ont., Canada L8S4KI. E-mail address:
[email protected] (M. Salahi).
0096-3003/$ - see front matter 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2006.05.092
M. Salahi, N. Mahdavi-Amiri / Applied Mathematics and Computation 183 (2006) 646–658
647
order algorithm [3] that has been considered in [10]. Before going into the details of the algorithm, first we briefly review the basic and unique results of IPMs. Throughout the paper we consider primal–dual IPMs to solve the following linear optimization (LO) problem: ðPÞ
minn fcT x : Ax ¼ b; x P 0g; x2R
where A 2 Rm·n satisfies rank(A) = m, b 2 Rm, c 2 Rn, and its dual problem ðDÞ
max fbT y : AT y þ s ¼ c; s P 0g:
y2Rm ;s2Rn
Without loss of generality [5] we may assume that both (P) and (D) satisfy the interior point condition (IPC), i.e., there exists an (x0, s0, y0) such that Ax0 ¼ b;
AT y 0 þ s0 ¼ c;
x0 > 0;
s0 > 0:
Under IPC, finding optimal solutions of (P) and (D) is equivalent to solving the following system: Ax ¼ b;
x P 0;
T
A y þ s ¼ c; Xs ¼ 0;
ð1Þ
s P 0;
where X = diag(x). The basic idea of primal–dual IPMs is to replace the third equation in (1) by the parameterized equation Xs = le, where e is the all one vector. This leads to the following system: Ax ¼ b; T
x P 0;
A y þ s ¼ c;
ð2Þ
s P 0;
Xs ¼ le: If the IPC holds, then for each l > 0, system (2) has a unique solution. This solution, denoted by (x(l), y(l), s(l)), is called the l-center of the primal–dual pair (P) and (D). The set of l-centers with all l > 0 gives the central path of (P) and (D) [2,8]. It has been shown that the limit of the central path (as l goes to zero) exists. Because the limit point satisfies the complementarity condition, it naturally yields optimal solutions for both (P) and (D), respectively [5]. One may consult [5,9] for more details on algorithmic developments. Now, we briefly describe the variation of Mehrotra’s second order predictor–corrector algorithm which is the focus of this paper. In the predictor step, the affine scaling search direction, A Dxa ¼ 0; AT Dy a þ Dsa ¼ 0; s Dxa þ x Dsa ¼ Xs
ð3Þ
is computed and the maximum step size in this direction is calculated so that ðx þ aa Dxa ; s þ aa Dsa Þ P 0: However, the algorithm does not take such a step right away. It uses the information from the predictor step and the second derivative of the primal dual trajectory to compute the centering direction given by A Dx ¼ 0; AT Dy þ Ds ¼ 0;
ð4Þ a
a
s Dx þ x Ds ¼ le Dx Ds ; where l is defined adaptively as 2 ga ga ; l¼ g n
648
M. Salahi, N. Mahdavi-Amiri / Applied Mathematics and Computation 183 (2006) 646–658
ga = (x + aa Dxa)T(s + aa Dsa) and g = xTs. Since (Dxa)T Dsa = 0, the previous relation implies 3
l ¼ ð1 aa Þ lg :
ð5Þ
Finally, the maximum step size a is computed so that the next iterate given by xðaÞ :¼ x þ a Dxa þ a2 Dx; yðaÞ :¼ y þ a Dy a þ a2 Dy; a
ð6Þ
2
sðaÞ :¼ s þ a Ds þ a Ds belong to a certain neighborhood of the central path. We call this step the corrector step of the algorithm. The polynomiality of this algorithm has been proved for monotone LCP in [10], not considering the need for adaptive update (needed for practical efficiency) of the centering parameter l at each iteration (they assume l = rklg at iteration k, where rmin 6 rk 6rmax, rmin and rmax are two constants in (0, 1)). They proved 3 that for the feasible case the algorithm has an O n2 log n worst case iteration complexity. We will use the numerical example presented in [7] and show that using the adaptive update of the barrier parameter can imply very small steps in order to keep the iterates in a certain neighborhood of the central path, and subsequently require a large number of iterations for the algorithm to stop with a desired solution. For practical purposes, we consider an adaptive update of barrier parameter, and to have control over a possible bad behavior of the algorithm, we introduce some safeguards. The rest of the paper is organized as follows. In Section 2, we describe the example and explain the possible bad behavior of the algorithm. Then in Section 3, we present the safeguard and establish the polynomial iteration complexity of the safeguarded algorithm. A slightly modified version of the safeguarded algorithm is presented in Section 4 which enjoys even a better iteration complexity. Computational results are given in Section 5. Finally, we close the paper with concluding remarks in Section 6. For self completeness we list three technical lemmas taken from [7,9] in Appendix A which are used frequently in the paper. For simplicity we use the following notation throughout the paper: I ¼ f1; . . . ; ng;
Iþ ¼ fi 2 IjDxai Dsai P 0g;
I ¼ fi 2 IjDxai Dsai < 0g;
F ¼ fðx; y; sÞ 2 Rn Rm Rn jAx ¼ b; AT y þ s ¼ c; x P 0; s P 0g; F0 ¼ fðx; y; sÞjAx ¼ b; AT y þ s ¼ c; x > 0; s > 0g; Dxc ¼ Dxa þ Dx;
Dy c ¼ Dy a þ Dy;
Dsc ¼ Dsa þ Ds;
lg ¼
xT s : n
2. Motivating example In this section first we introduce the neighborhood of the central path that will be used throughout the paper. Then, we describe a numerical example showing that variation of Mehrotra’s second order predictor–corrector algorithm as described in the previous section might take very small steps to keep the iterates in the defined neighborhood of the central path, and hence take many iterations to convergence. The example indicates that Mehrotra’s adaptive updating scheme of the centering (barrier) parameter and correction terms in system (4) need to be combined with certain safeguards to get a warranted step size at each iteration and achieve a polynomial iteration complexity. Knowing that most efficient IPMs based solvers work in the negative infinity norm neighborhood defined by 0 N 1 ðcÞ :¼ fðx; y; sÞ 2 F : xi si P clg 8i 2 Ig;
where c 2 (0, 1) is a constant independent of n, we consider algorithms that are working in N 1 ðcÞ.
ð7Þ
M. Salahi, N. Mahdavi-Amiri / Applied Mathematics and Computation 183 (2006) 646–658
649
Let us consider the following simple LO: min s:t:
x2 0 6 x1 6 1;
ð8Þ
0 6 x2 6 1 þ dx1 : Let the algorithm start from the following feasible points in the neighborhood N 1 ð0:098Þ: x0 ¼ ð0:03; 0:9; 0:97; 0:1022Þ;
s0 ¼ ð5:552; 1; 5:7; 2Þ;
y 0 ¼ ð5:7; 2Þ;
where d = 0.074. For the given starting point, by using identical step sizes for both primal and dual problems,1 in the first iteration the maximum step size in the predictor step is aa = 0.95 while the maximum step size in the corrector step is Oð102 Þ and even getting smaller for later iterations. To explain our observation, we examine the constraints xi ðaÞsi ðaÞ P clg ðaÞ
8i 2 I
ð9Þ
for c = 0.098 that keep the next iterate in N 1 ð0:098Þ neighborhood. By expanding inequality (9) and reordering we obtain: ð1 aÞxi si þ a2 ð1 cÞl þ a3 ðDxai Dsi þ Dx Dsai Þ þ a4 Dxi Dsi P cð1 aÞlg
8i 2 I:
Note that for the given staring point, x1s1 0.098lg is zero, which means that the first coordinate of xs is on the boundary of the neighborhood. We also have Dxa1 Ds1 þ Dx1 Dsa1 < 0, Dx1 Ds1 < 0, l ¼ Oð104 Þ, and the absolute value of Dxa1 Ds1 þ Dx1 Dsa1 dominates l. Incorporating all these information into (9) implies that the algorithm requires a very small step ðOð102 ÞÞ to satisfy (9). This phenomenon may have been caused by: • an aggressive update of centering parameter l using (5); • the usage of the correction terms in the centering system of equations. To alleviate the aforementioned phenomenon, we propose the following remedies: • Cutting the maximum step size in the predictor step if it is greater than a certain threshold. This might prevent the algorithm from an aggressive update. • Removing the correction terms in the corrector system of equations. For this specific example, these ideas help us to solve the difficulty that might arise. However, in general, removing the correction terms is not as effective as cutting the maximum step size in the predictor step and using (5) to compute l. These observations motivate us to introduce some safeguards that will help us to have control over the minimal warranted step size both theoretically and practically. In our safeguard we simply cut the maximum step size in the predictor step if it is above a certain threshold. 3. Safeguarded algorithm In this section we give a condition on the maximum step size in the predictor step, aa, violation of which might imply a very small step size for the corrector step. This motivates us to cut aa when it violates the given condition. Motivated by this observation we will modify the algorithm as presented at the end of this section. The following lemma gives a lower bound for the maximum step size in the predictor step and the proof can be found in [6,7]. This result indicates that there always exists a positive step size in the predictor step. a a a Theorem 3.1. Suppose that the current iterate ðx; y; sÞ 2 N 1 ðcÞ and (Dx , Dy , Ds ) be the solution of (3). Then the maximum feasible step size, aa, such that (x(aa), s(aa)) P 0, satisfies
1
In practice one might use different step sizes in primal and dual spaces. By slightly modifying the starting point of the example, one can realize similar results when taking different step sizes.
650
M. Salahi, N. Mahdavi-Amiri / Applied Mathematics and Computation 183 (2006) 646–658
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c2 þ cn 2c : aa P n 2
ð10Þ
The following technical lemmas are useful in the step size estimation for the corrector step of the new algorithm. a a a Lemma 3.2. Suppose that the current iterate ðx; y; sÞ 2 N 1 ðcÞ and (Dx , Dy , Ds ) be the solution of (3) and (Dx, Dy, Ds) be the solution of (4) with l P 0. Then,
kDxa Dsk 6
!12 2 1 l 1 1 l 3 þ þ n2 lg þ 2c lg 32 16c 4clg
kDx Dsa k 6
!12 2 1 l 1 1 l 3 þ þ þ n2 lg : 2c lg 32 16c 4clg
and
1
Proof. Let D ¼ ðX 1 SÞ2 . Then, by multiplying the third equation of (3) by ðXSÞ
12
we have
1
D Dxa þ D1 Dsa ¼ ðXsÞ2 : Therefore, by squaring both sides of this equation and due to the orthogonality of Dxa and Dsa one has pffiffiffiffiffiffiffi pffiffiffiffiffiffiffi ð11Þ kD Dxa k 6 nlg and kD Dsa k 6 nlg : By doing the same procedure for the third equation of (4) one has 1
D Dx þ D1 Ds ¼ ðXsÞ 2 ðl Dxa Dsa Þ: Now, using Lemmas A.1 and A.2 and due to the orthogonality of Dx and Ds we obtain: !12 2 1 2 nlg n2 lg nl 2 nl 1 l 1 1 l pffiffiffiffiffi kD Dxk 6 þ þ þ þ 6 þ þ n lg 2c lg 32 16c 4clg clg 16 16c 2c
ð12Þ
and kD
1
2 1 nlg n2 lg nl 2 nl þ þ Dsk 6 þ 6 clg 16 16c 2c
!12 2 1 l 1 1 l pffiffiffiffiffi þ þ þ n lg : 2c lg 32 16c 4clg
ð13Þ
Finally, using the fact that D is diagonal and inequalities (11)–(13) we have !12 2 1 l 1 1 l 3 þ kDxa Dsk 6 kD Dxa kkD1 Dsk 6 þ þ n2 lg 2c lg 32 16c 4clg and analogously one has the second statement of the lemma, which completes the proof of the lemma.
h
In the following lemma we give an upper bound for kDx Dsk to be used in our later analysis. One can find the details of proof in [6,7]. Lemma 3.3. Suppose that the current iterate is ðx; y; sÞ 2 N 1 ðcÞ and let (Dx, Dy, Ds) be the solution of (4), where l P 0. Then we have 2 nlg n2 lg nl 3 nl 2 þ þ kDx Dsk 6 2 þ : clg 16 16c 2c 1
Proof. If we multiply the third equation of (4) by ðXSÞ 2 , then by Lemma A.3 we have the desired result.
h
M. Salahi, N. Mahdavi-Amiri / Applied Mathematics and Computation 183 (2006) 646–658
651
Now, in the following theorem we give a condition on aa violation of which might imply a very small step size in the corrector step of the algorithm. a a a Theorem 3.4. Suppose that the current iterate ðx; y; sÞ 2 N 1 ðcÞ and (Dx , Dy , Ds ) be the solution of (3) and (Dx, Dy, Ds) be the solution of (4) with l as defined by (5). Then, for aa satisfying
ct aa < 1 ð1 cÞ where
t ¼ max i2Iþ
13 ð14Þ
;
Dxai Dsai ; x i si
ð15Þ
the algorithm always takes a step with positive step size in the corrector step. Proof. The goal is to determine maximum a 2 (0, 1] such that (9) holds. To do so, one has xi ðaÞsi ðaÞ ¼ ðxi þ a Dxai þ a2 Dxi Þðsi þ a Dsai þ a2 Dsi Þ ¼ ð1 aÞxi si þ a2 l a2 Dxai Dsai þ a2 Dxai Dsai þ a3 ðDxai Dsi þ Dx Dsai Þ þ a4 Dxi Dsi ¼ ð1 aÞxi si þ a2 l þ a3 ðDxai Dsi þ Dx Dsai Þ þ a4 Dxi Dsi :
ð16Þ
We also have lg ðaÞ ¼ ð1 aÞlg þ a2 l: Therefore, (9) is equivalent to
ð1 aÞxi si þ a2 l þ a3 Dxai Dsi þ Dx Dsai þ a4 Dxi Dsi P cð1 aÞlg þ ca2 l
ð17Þ
8i 2 I:
The worst case for this inequality happens when xi si = clg, Dxai Dsi þ Dx Dsai < 0 and Dxi Dsi < 0. Then, the previous inequality is equivalent to
ð1 cÞl þ a Dxai Dsi þ Dx Dsai þ a2 Dxi Dsi P 0: Violating this inequality as it happened in the example (8) of the previous section might result in a very small step size in the corrector step. On the other hand, one can notice that this inequality involves directions of the centering step. In what follows we derive a slightly weaker condition, not involving the centering direction, which can be verified after the predictor step and this way one saves a backsolve which needs to be done in case of violating the previous condition. Since (Dxa)T Dsa = 0, then Iþ 6¼ ;. Now, for i 2 Iþ we can write relation (16) as xi ðaÞsi ðaÞ ¼ ð1 aÞxi si þ a2 l a2 Dxai Dsai þ a2 Dxai Dsai þ a3 ðDxai Dsi þ Dxi Dsai Þ þ a4 Dxi Dsi P ð1 aÞxi si þ a2 l a2 Dxai Dsai þ a3 ðDxai Dsai þ Dxai Dsi þ Dxi Dsai þ Dxi Dsi Þ ¼ ð1 aÞxi si þ a2 l a2 Dxai Dsai þ a3 Dxci Dsci ; where the inequality follows from the assumption we made for the worst case analysis, namely Dxi Dsi < 0. Finally, for i 2 Iþ , (9) holds whenever ð1 aÞxi si þ a2 l a2 Dxai Dsai þ a3 Dxci Dsci P cð1 aÞlg þ ca2 l: Since Dxai Dsai 6 txi si by the definition of t given in (15), it is sufficient to have a P 0 for which ð1 a a2 tÞxi si þ a2 l þ a3 Dxci Dsci P cð1 aÞlg þ ca2 l: Using the fact that xisi P clg and t 6 14 by Lemma A.1, the previous inequality holds for a 6 1þ2pffiffi2 whenever ð1 a a2 tÞclg þ a2 l þ a3 Dxci Dsci P cð1 aÞlg þ ca2 l
652
M. Salahi, N. Mahdavi-Amiri / Applied Mathematics and Computation 183 (2006) 646–658
or ð1 cÞl ctlg þ a Dxci Dsci P 0: We also assume that in the worst case Dxci Dsci < 0, then one has to have (1 c)(1 aa)3 ct > 0. This def 13 ct . The corresponding inequalities in (9) for i 2 I also holds for these initely holds whenever aa < 1 1c h
values of aa, which completes the proof.
Since by Lemma A.1 t varies in ½0 14, then in order to guarantee a unique lower bound for the maximum step size in thecorrector step for the purpose of polynomial iteration complexity of the algorithm, we let 1 3 c aa ¼ 1 2ð1cÞ whenever the maximum step size in the corrector step is below certain threshold. In the following lemma we give the lower for the maximum step size in the corrector step for this specific choice. bound 1 Note also that for aa ¼ 1
c 2ð1cÞ
3
c by using (5) one has l ¼ 2ð1cÞ lg . The following corollaries which follow
from Lemmas 3.2 and 3.3 are useful in the next theorem. c Corollary 3.5. Let l ¼ 2ð1cÞ lg . Then 3
n2 kDx Dsk 6 pffiffiffi lg 2 c a
and 3
n2 kDx Ds k 6 pffiffiffi lg : 2 c a
c Corollary 3.6. Let l ¼ 2ð1cÞ lg . Then
n2 kDx Dsk 6 pffiffiffi lg : 8 2c a a a Theorem 3.7. Suppose that the current iterate ðx; y; sÞ 2 N 1 ðcÞ and (Dx , Dy , Ds ) be the solution of (3) and c (Dx, Dy, Ds) be the solution of (4) with l ¼ 2ð1cÞ lg . Then, 3
aP
c2 3
3n2
:
Proof. The goal is to determine maximum step size a 2 (0, 1] in the corrector step such that (9) holds. Following the similar analysis of the previous theorem it is sufficient to have ð1 cÞl þ aðDxai Dsi þ Dxi Dsai Þ þ a2 Dxi Dsi P 0: Using Corollaries 3.5 and 3.6 it is sufficient to have 1 1 n2 3 clg pffiffiffi an2 lg pffiffiffi a2 P 0: 2 c 8 2c 3
This inequality definitely holds for a ¼
c2 3
, which completes the proof.
h
3n2
Remark 3.8. When the algorithm uses the adaptive update of barrier parameter, a small aa might imply a marginal reduction of the duality gap given by (17). Therefore, for aa small, for example aa < 0.1, the new algorithm uses our new safeguard. Based on the aforementioned discussions, we can now outline our safeguard based second order Mehrotratype predictor–corrector algorithm as Algorithm 1.
M. Salahi, N. Mahdavi-Amiri / Applied Mathematics and Computation 183 (2006) 646–658
653
Algorithm 1 (Safeguarded algorithm). Input: A proximity parameters c 2 ð0; 13Þ; an accuracy parameter > 0; ðx0 ; y 0 ; s0 Þ 2 N 1 ðcÞ. begin while xTs P do begin Predictor Step Solve (3) and compute the maximum step size aa such that ðxðaa Þ; yðaa Þ; sðaa ÞÞ 2 F; end begin Corrector step If aa P 0.1, then solve (4) with l = (1 aa)3lg and compute the maximum step size a such that (x(a), y(a), s(a)) given by (6) belongs to N 1 ðcÞ; 3
If a <
c2 3
3n2
c , then solve (4) with l ¼ 2ð1cÞ lg and compute the maximum step size a such that
ðxðaÞ; yðaÞ; sðaÞÞ 2 N 1 ðcÞ; end else c Solve (4) with l ¼ 2ð1cÞ lg and compute the maximum step size a such that ðxðaÞ; yðaÞ; sðaÞÞ 2 N 1 ðcÞ; end Set (x, y, s) = (x(a), y(a), s(a)). end end In the following theorem we give an upper bound for the maximum number of iterations that is required by Algorithm 1 to find an -approximate solution. Theorem 3.9. Algorithm 1 stops after at most ! ðx0 ÞT s0 3 O n2 log number of iterations with a solution for which xTs 6 . 3
Proof. If aa P 0.1 and a P
c2
, then using (17) we have ! 3 2c2 2 l lg ðaÞ ¼ 1 a þ a l 6 1 lg : 3 lg g 25n2 3
3n2
3
If aa P 0.1 and a < lg ðaÞ ¼
c2 3
, then
3n2
l 1aþa l 6 lg g 2
3
1
c2 ð2 3cÞ 3
6ð1 cÞn2
! lg :
Finally, if aa < 0.1, then again one has ! 3 c2 ð2 3cÞ 2 l lg ðaÞ ¼ 1 a þ a l 6 1 lg ; 3 lg g 6ð1 cÞn2 which completes the proof of the theorem. h
654
M. Salahi, N. Mahdavi-Amiri / Applied Mathematics and Computation 183 (2006) 646–658
Remark 3.10. It is worth mentioning that Algorithm 1 has a better iteration complexity than its first order analog in [6,7]. 4. A modification of Algorithm 1 In this section we present a slight modification of Algorithm 1 that reduces the iteration complexity significantly. The motivation for this modification comes from the following lemma that reduces the bound on lemmas significantly and leads to a slight modification of (6). Lemma 4.1. For i 2 I and for any 0 < aa 6 1 we have 1 1 Dxai Dsai 6 1 xi si : aa aa
ð18Þ
Proof. For the maximum feasible step size in the predictor step, aa, we have xi ðaa Þsi ðaa Þ P 0;
i ¼ 1; . . . ; n:
This is equivalent to ð1 aa Þxi si þ a2a Dxai Dsai P 0;
i ¼ 1; . . . ; n:
For i 2 I the statement of the lemma follows. h Based on the aforementioned observation we modify the centering system of equations by introducing a factor of aa in the right hand side of the third equation of the centering step. The new system becomes A Dx ¼ 0; AT Dy þ Ds ¼ 0; s Dx þ x Ds ¼ le aa Dxa Dsa ;
ð19Þ
where the centering parameter l is defined by (5). In other words, in the new algorithm we solve system (19) instead of system (4) in the centering step. We also change the corrector step given in (6) to the following: pffiffiffiffiffi xðaÞ :¼ x þ a aa Dxa þ a2 Dx; pffiffiffiffiffi ð20Þ yðaÞ :¼ y þ a aa Dy a þ a2 Dy; pffiffiffiffiffi a 2 sðaÞ :¼ s þ a aa Ds þ a Ds: Analogous to Lemma 3.3 we have the following bound for kDx Dsk, which is stronger than the bound given in Lemma 3.3. Lemma 4.2 [7]. Suppose that the current iterate ðx; y; sÞ 2 N 1 ðcÞ, 0 < aa 6 1 and (Dx, Dy, Ds) be the solution of (19). Then we have ! 2 2 3 1 l a 4 4a a l a þ aa kDx Dsk 6 2 2 þ þ nlg : c lg 2c lg 16 a a a Lemma 4.3. Suppose that the current iterate ðx; y; sÞ 2 N 1 ðcÞ and (Dx , Dy , Ds ) be the solution of (3) and (Dx, Dy, Ds) be the solution of (4) with l as defined by (5). Then,
kDxa Dsk 6
!12 2 1 l aa l 4 4aa þ a2a þ þ nlg c lg 2c lg 16
M. Salahi, N. Mahdavi-Amiri / Applied Mathematics and Computation 183 (2006) 646–658
655
and kDx Dsa k 6
!12 2 1 l aa l 4 4aa þ a2a þ þ nlg : c lg 2c lg 16
Analogous to the previous section, in the following theorem we give a condition on aa violation of which might imply a very small step size in the corrector step. The proof is similar to the proof of Theorem 3.4. a a a Theorem 4.4. Suppose that the current iterate ðx; y; sÞ 2 N 1 ðcÞ and (Dx , Dy , Ds ) be the solution of (3) and (Dx, Dy, Ds) be the solution of (19) with l as defined by (5). Then, for aa satisfying 1 caa t 3 aa < 1 ð21Þ ð1 cÞ
the modified algorithm always takes a step with positive step size in the corrector step. Following the 13 discussion of the previous section, for the sake of polynomial complexity we let c c aa ¼ 1 2ð1cÞ when ac is below certain threshold. This further implies l ¼ 2ð1cÞ lg using (5). c Corollary 4.5. Let l ¼ 2ð1cÞ lg . Then pffiffiffi 7n kDxa Dsk 6 pffiffiffiffiffi lg 2 2c and pffiffiffi 7n a kDx Ds k 6 pffiffiffiffiffi lg : 2 2c c lg . Then Corollary 4.6. Let l ¼ 2ð1cÞ 7n kDx Dsk 6 pffiffiffi lg : 16 2c
a a a Lemma 4.7. Suppose that the current iterate ðx; y; sÞ 2 N 1 ðcÞ and (Dx , Dy , Ds ) be the solution of (3) and c (Dx, Dy, Ds) be the solution of (19) with l ¼ 2ð1cÞ lg . Then, 3
c2 a P pffiffiffiffiffi : 4 aa n Proof. The goal is to determine maximum a 2 (0, 1] such that (9) holds. Following the similar analysis of Theorem it is sufficient to have
pffiffiffiffiffi ð1 cÞl þ a aa Dxai Dsi þ Dx Dsai þ a2 Dxi Dsi P 0: Using Corollaries 4.5 and 4.6 it is sufficient to have pffiffiffi 1 7 7 pffiffiffiffiffi cl pffiffiffiffiffi aa anlg pffiffiffi a2 n P 0: 2 g 2c 16 2c 3 2
This inequality definitely holds for a ¼ 4pcffiffiffi , which completes the proof. aa n Theorem 4.8. The modified version of Algorithm 1 stops after at most ! ðx0 ÞT s0 O n log number of iterations with a solution for which xTs 6 .
h
656
M. Salahi, N. Mahdavi-Amiri / Applied Mathematics and Computation 183 (2006) 646–658
Proof. The proof is analogous to the proof of Theorem 3.9.
h
5. Computational results In this section we report some limited encouraging computational results for the test problems given in Table 1 that are taken from the standard NETLIB test repository [4]. The implementation of our new safeguarded algorithms is carried out by modifying some subroutines of the McIPM software package [12], a self-dual embedding implementation of IPMs for LO, on Matlab 6.5, using a Pentium 4, 2.53 GHz, and 512 MB ram. We let c = 0.001 in our implementation and results are summarized in Table 2. We have implemented a combination of Algorithm 1 and its modified version in the following way. If the affine scaling is doing well, aa P 0.1, then we use system (4), otherwise we use system (19). The following abbreviations are used for different variants of McIPM software package: • PMcIPM denotes a version of McIPM which is using Mehrotra’s updating scheme given by (5) combined with certain heuristics. • SMcIPM denotes a version of McIPM which is using the safeguarded algorithm. Computational results show that our new safeguarded algorithm not only enjoys a polynomial iteration complexity but also requires less number of iteration for most of the problems. The difference sometimes can be dramatic, for example for problems ‘degen3’, ‘pds20’ and ‘pilot’. The running times of the test problems are not reported, since significant differences occur only if the number of iterations are significantly different, for example for ‘degen3’ and ‘pds-20’. Table 1 Test problems Problem
Rows
Columns
Problem
Rows
Columns
cre-a cre-b cre-c cre-d degen3 ken-07 ken-11 ken-13 ken-18 osa-07 osa-14
3517 9649 3069 8927 1503 2427 14,695 28,633 105,128 1119 2338
4067 72,477 3678 69,980 2604 3602 21,349 42,659 154,699 23,949 52,460
osa-30 osa-60 perold pds-02 pds-06 pds-10 pds-20 pilot scsd1 scsd6 scsd8
4351 10,281 625 2954 9882 16,559 33,875 1441 760 1350 3750
100,024 232,966 1530 7535 28,655 48,763 105,728 4657 3648 5666 11,334
Table 2 Comparison of iterations number for variants of McIPM on the test problems Problem
PMMcIPM
SMcIPM
Problem
PMMcIPM
SMcIPM
cre-a cre-b cre-c cre-d degen3 ken-07 ken-11 ken-13 ken-18 osa-07 osa-14
27 34 31 33 22 18 20 29 37 31 39
27 34 30 33 14 17 20 28 36 30 38
osa-30 osa-60 perold pds-02 pds-06 pds-10 pds-20 pilot scsd1 scsd6 scsd8
44 51 45 34 52 73 96 61 12 13 12
42 48 43 33 50 74 85 51 11 12 11
M. Salahi, N. Mahdavi-Amiri / Applied Mathematics and Computation 183 (2006) 646–658
657
6. Conclusions In this paper a variation of Mehrotra’s second order predictor–corrector algorithm is analyzed. By a numerical example it is shown that the algorithm might suffer from some drawbacks which motivated us to propose some safeguards. We proved that the safeguarded algorithm has polynomial iteration complexity. It is also proved that a second modification of the algorithm enjoys even a better iteration complexity compared to the similar algorithm in [10]. Finally, some limited encouraging computational results using the McIPM software package is reported. Acknowledgments The research of the first author is supported by the NSERC Discovery Grant and by MITACS through Prof. Tama´s Terlaky. Appendix A In this section we present some technical lemmas that are used frequently throughout the paper. These lemmas are quoted from [7]. Lemma A.1. Let (Dxa, Dya, Dsa) be the solution of (3), then Dxai Dsai 6
xi si 4
8i 2 Iþ :
Proof. By Eq. (3) for i 2 Iþ we have si Dxai þ xi Dsai ¼ xi si : If we divide this equation by xisi we get Dxai Dsai þ ¼ 1: xi si Since Dxai Dsai > 0, this equality implies that both Dxai < 0 and Dsai < 0. Then one can easily prove that Dxai Dsai 6
xi si : 4
Lemma A.2. Let (Dxa, Dya, Dsa) be the solution of (3), then we have X
Dxai Dsai ¼
i2Iþ
X
jDxai Dsai j 6
i2I
xT s : 4
Proof. The proof is a direct consequence of the previous lemma.
h
The following technical lemma is recalled several times throughout the paper and the proof can be found in [9]. Lemma A.3. Let u and v be any two vectors in Rn with uTv P 0. Then 3
2
kUVek 6 2 2 ku þ vk ; where U = diag(u) and V = diag(v).
658
M. Salahi, N. Mahdavi-Amiri / Applied Mathematics and Computation 183 (2006) 646–658
References [1] J. Czyzyk, S. Mehrtotra, M. Wagner, S.J. Wright, PCx: An interior-point code for linear programming, Optimization Methods and Software 11–12 (1999) 397–430. [2] N. Megiddo, Pathways to the optimal set in linear programming, in: N. Megiddo (Ed.), Progress in Mathematical Programming: Interior Point and Related Methods, Springer-Verlag, New York, 1986, pp. 131–158 (Identical version in: Proceedings of the 6th Mathematical Programming Symposium of Japan, Nagoya, Japan, 1–35). [3] S. Mehrotra, On the implementation of a (primal–dual) interior point method, SIAM Journal on Optimization 2 (1992) 575–601. [4] NETLIB: A repository of linear programming test problems. Available from: . [5] C. Roos, T. Terlaky, J.-Ph. Vial, Theory and Algorithms for Linear Optimization: An Interior Point Approach, John Wiley & Sons, Chichester, UK, 1997. [6] M. Salahi, New Adaptive Interior Point Algorithms for Linear Optimization, Ph.D. Thesis, Advanced Optimization Lab., McMaster University, Canada, March 2006. [7] M. Salahi, J. Peng, T. Terlaky, On Mehrtora type predictor–corrector algorithms, Technical Report 2005/4, Advanced Optimization Lab, Department of Computing and Software, McMaster University, http://www.cas.mcmaster.ca/~oplab/publication, SIAM Journal on Optimization, submitted for publication. [8] G. Sonnevend, An ‘‘analytic center’’ for polyhedrons and new classes of global algorithms for linear (smooth, convex) programming, in: A. Pre´kopa, J. Szelezsa´n, B. Strazicky (Eds.), System Modeling and Optimization: Proceeding of the 12th IFIP Conference, Budapest, Hungary, September 1985, Lecture Notes in Control and Information Sciences, vol. 84, Springer-Verlag, Berlin, 1986, pp. 866–876. [9] S.J. Wright, Primal-Dual Interior-Point Methods, SIAM, Philadelphia, USA, 1997. [10] Y. Zhang, D. Zhang, On the polynomiality of the Mehrotra-type predictor–corrector interior point algorithms, Mathematical Programming 68 (1995) 303–317. [11] Y. Zhang, Solving large-scale linear programms by interior point methods under the Matlab enviroment, Optimization Methods and Software 10 (1999) 1–31. [12] X. Zhu, J. Peng, T. Terlaky, G. Zhang, On implementing self-regular proximity based feasible IPMs, Technical Report 2003/2, Advanced Optimization Lab, Department of Computing and Software, McMaster University, http://www.cas.mcmaster.ca/~oplab/ publication, INFORMS Journal on Computing, submitted for publication.
