# A new logarithmic-quadratic proximal method for nonlinear complementarity problems

## A new logarithmic-quadratic proximal method for nonlinear complementarity problems

Applied Mathematics and Computation 215 (2009) 695–706 Contents lists available at ScienceDirect Applied Mathematics and Computation journal homepag...

Applied Mathematics and Computation 215 (2009) 695–706

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

A new logarithmic-quadratic proximal method for nonlinear complementarity problems Abdellah Bnouhachem a,b, Muhammad Aslam Noor c,*, Mohamed Khalfaoui d, Sheng Zhaohan a a

School of Management Science and Engineering, Nanjing University, Nanjing 210093, PR China Ibn Zohr University, ENSA, BP 32/S, Agadir, Morocco Mathematics Department, COMSATS Institute of Information Technology, Sector-H-8/1, Islamabad 44000, Pakistan d Ibn Zohr University, Présidence, BP 32/S, Agadir, Morocco b c

a r t i c l e

i n f o

a b s t r a c t In this paper, we propose a new modiﬁed logarithmic-quadratic proximal (LQP) method for solving nonlinear complementarity problems (NCP). We suggest using a prediction-correction method to solve NCP. The predictor is obtained via solving the LQP system approximately under signiﬁcantly relaxed accuracy criterion and the new iterate is computed by using a new step size ak . Under suitable conditions, we prove that the new method is globally convergent. We report preliminary computational results to illustrate the efﬁciency of the proposed method. This new method can be considered as a signiﬁcant reﬁnement of the previously known methods for solving nonlinear complementarity problems. Ó 2009 Elsevier Inc. All rights reserved.

Keywords: Nonlinear complementarity problems Pseudomonotone operators Interior proximal methods

1. Introduction The nonlinear complementarity problem (NCP) is to determine a vector x 2 Rn such that

x P 0;

FðxÞ P 0 and xT FðxÞ ¼ 0;

ð1:1Þ

n

where F is a nonlinear mapping from R into itself. For the applications, formulation and numerical method of the complementarity problems and related variational inequalities, see Refs. [2–29] and the references therein. It is well known that NCP can be alternatively formulated as ﬁnding the zero of an appropriately maximal monotone operator T ¼ F þ N Rnþ , namely, ﬁnd x 2 Rnþ such that 0 2 Tðx Þ, where N Rnþ ð:Þ is the normal cone operator to Rnþ deﬁned by

NRnþ ðxÞ :¼

( y : yT ðv  xÞ 6 0; ;;



8v 2 Rnþ ; if

x 2 Rnþ ;

otherwise:

A well known method to ﬁnd the zero of a maximal monotone operator T is the proximal point algorithm (PPA). It was ﬁrst introduced by Martinet  and further extended by Rockafellar  who showed how the PPA could be applied to convex programs, convex saddle point problems and variational inequality problems. PPA starts with any vector x0 2 Rnþ and bk P b > 0, iteratively updates xkþ1 conforming the following problem:

0 2 bk TðxÞ þ rx qðx; xk Þ;

ð1:2Þ

where

qðx; xk Þ ¼

1 kx  xk k2 2

* Corresponding author. E-mail addresses: [email protected] (A. Bnouhachem), [email protected], [email protected] (M.A. Noor). 0096-3003/\$ - see front matter Ó 2009 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2009.05.042

ð1:3Þ

696

A. Bnouhachem et al. / Applied Mathematics and Computation 215 (2009) 695–706

is a quadratic function of x. Motivation for studying the algorithms of problem (1.2) could be found in several studies [13,23,27], in place of usual quadratic term where many researchers have used some nonlinear functions rðx; xk Þ. For instance, we quote reference  for the iterative schemes of the form (1.2) using the Bregman-based functional instead of (1.3). The inexact version of the proximal point algorithm (PPA)  generates iteratively sequence fxk g  Rnþ satisfying

ek 2 ðxkþ1  xk Þ þ bk Tðxkþ1 Þ

ð1:4Þ

where ek 2 Rn is the error term. In the last decade, there have been increasing interests in constructing new error criteria of PPA [7,9,14,23] and in using growth condition to ensure convergence behavior the proposed iterative algorithms. Rockafellar  proposed the ﬁrst relative error criterion of PPA that all the relaxation factors are summable. Eckstein  suppose that 1 X

kek k < þ1 and

k¼1

1 X hek ; xk i < þ1:

ð1:5Þ

k¼1

On the other hand, Han and He  proved the convergence under the following accuracy criterion

kek k 6 gk kxk  xk1 k with

þ1 X

g2k < þ1:

ð1:6Þ

k¼0

Furthermore, Da Silva et al.  and Solodov and Svaiter [25,26] proposed some new accurate criteria for PPA. Their criteria require only that supkP0 gk < 1. Recently, Auslender et al.  have proposed a new type of proximal interior method through replacing the quadratic function (1.3) by d/ ðx; xk Þ which could be deﬁned as

d/ ðx; yÞ ¼

n X

y2j /ðy1 j xj Þ:

j¼1

The fundamental difference here is that the term d/ is used to force the iterates fxkþ1 g to stay in the interior of the nonnegative orthant Rnþþ . Among the possible choices of /, there exists a particular one which enjoys several attractive properties for developing efﬁcient algorithms to solve NCP. Let m > l > 0 be given ﬁxed parameters, and deﬁne

(

m ðt

2

/ðtÞ ¼

 1Þ2 þ lðt  log t  1Þ if

þ1

t>0

otherwise:

In , Auslender et al. used a very special logarithmic-quadratic proximal (LQP) method (with m ¼ 2; l ¼ 1) for solving variational inequalities over polyhedra. In order to ensure convergence, they have suggested to use the accuracy criterion of type (1.5). In this paper, we considered m ¼ 1 and l 2 ð0; 1Þ, so the function / is deﬁned by

( /ðtÞ ¼

1 ðt 2

 1Þ2 þ lðt  log t  1Þ if

þ1

t>0

otherwise:

Then the problem (1.2) becomes for given xk 2 Rnþ and bk P b > 0, the new iterate xkþ1 is unique solution of the following setvalued equation:

0 2 bk TðxÞ þ rx Q ðx; xk Þ;

ð1:7Þ

8    2  k n 2 > < 1 kx  xk k2 þ l P xk log xj þ xj xk  xk if x 2 Rnþþ j j j xj Q ðx; xk Þ ¼ 2 j¼1 > : þ1 otherwise:

ð1:8Þ

where

It is easy to see that

rx Q ðx; xk Þ ¼ x  ð1  lÞxk  lX 2k x1 ;  where X k ¼ diag xk1 ; xk2 ; . . . ; xkn , x1 is an n-vector whose jth element is 1=xj and l 2 ð0; 1Þ is a given constant. Then the problem (1.7) and (1.8) is equivalent to the following systems of nonlinear equations

bk FðxÞ þ x  ð1  lÞxk  lX 2k x1 ¼ 0:

ð1:9Þ

In Auslender, Teboulle and Ben-Tiba’s original work , it is required that the subproblem should be solved exactly. To release the difﬁculty in Auslender’s problem, very recently, He et al. , Bnouhachem , Noor and Bnouhachem ,

A. Bnouhachem et al. / Applied Mathematics and Computation 215 (2009) 695–706

697

and Xu et al.  introduced some LQP based prediction-correction methods which do not suffer from the above difﬁculty and make the LQP method more practical. Each iteration of the above methods contains a prediction and a correction, the predictor is obtained via solving the LQP system approximately under signiﬁcantly relaxed accuracy criterion and the new iterate is computed directly by an explicit formula derived from the original LQP method for , while the new iterate is computed by using the projection operator for [5,21,28]. Inspired and motivated by the above research, we suggest and analyze a new LQP method for solving nonlinear complementarity problems (1.1) by using a new step size ak which provides a signiﬁcant reﬁnement and improvement of the method of . We also study the global convergence of the proposed modiﬁed LQP method under some mild conditions. Some preliminary computational results are given to illustrate the efﬁciency of the proposed method. Throughout this paper we assume that F is continuously differentiable and pseudomonotone with respect to Rþ n and the solution set of (1.1), denoted by X , is nonempty.

2. Preliminaries We ﬁrst recall some important results which are main tool for our theoretical analysis. We denote P Rnþ ð:Þ as the projection under the Euclidean norm, i.e.,

  PRnþ ðzÞ ¼ min kz  xkj x 2 Rnþ : A basic property of the mapping of projection is

hy  PRnþ ðyÞ; PRnþ ðyÞ  xi P 0;

8y 2 Rn ;

8x 2 Rnþ :

ð2:1Þ

From (2.1), it is easy to verify that

kPRnþ ðv Þ  uk2 6 kv  uk2  kv  PRnþ ðv Þk2 ;

8v 2 Rn ; u 2 Rnþ :

ð2:2Þ

Deﬁnition 2.1. The operator F : Rn ! Rn is said to be pseudomonotone, if

8u; v 2 Rn ;

hv  u; FðuÞi P 0 ) hv  u; Fðv Þi P 0:

The following lemma is similar to Lemma 2 in . For the sake of completeness and to convey an idea of the technique involved, we include its proof. Lemma 2.1. For given xk > 0 and q 2 Rn , let x be the positive solution of the following equation:

q þ x  ð1  lÞxk  lX 2k x1 ¼ 0;

ð2:3Þ

then for any y P 0 we have

hy  x; qi P

 1l 1þl kx  yk2  kxk  yk2 þ kxk  xk2 : 2 2

ð2:4Þ

Proof. Since x > 0, xk > 0 and y P 0, we have

 2 yi xki =xi P yi ð2xki  xi Þ;

i ¼ 1; . . . ; n:

ð2:5Þ

By a manipulation, it follows from (2.3) that for i ¼ 1; . . . ; n,

n  2 o  2  ðyi  xi Þqi ¼ ðxi  yi Þ xi  ð1  lÞxki  l xki =xi P ðxi Þ2  ð1  lÞxi xki  l xki  xi yi þ ð1  lÞxki yi þ lyi 2xki  xi  2  2  1  l  k 2 1þl þ ðxi  yi Þ2  xki  yi xi  x i : ¼ ðxi Þ2  ð1  lÞxi xki  l xki  ð1 þ lÞxi yi þ ð1 þ lÞxki yi ¼ 2 2 Hence (2.4) holds and the proof is completed.

h

  xk 2 Rnþ , bk > 0, and xkþ1 ðaÞ is deﬁned by xkþ1 ðak Þ ¼ sxk þ ð1  sÞPRnþ hLemma 2.2. iLet x 2 X , 0 < s < 1; 0 < l < 1, ~ k k bk ~ xk  a1þ Fð x Þ , then we have l

2ak bk k 2ak bk k kxkþ1 ðak Þ  x k2 6 kxk  x k2  ð1  sÞ kxk  xk ðak Þk2 þ hx ðak Þ  xk ; Fð~xk Þi þ hx  x ; Fð~xk Þi ; 1þl 1þl

h i k bk ~k where xk ðak Þ ¼ PRnþ xk  a1þ l Fðx Þ .

698

A. Bnouhachem et al. / Applied Mathematics and Computation 215 (2009) 695–706

h i k bk ~k Proof. Since x 2 X  Rnþ and xk ðak Þ ¼ PRnþ xk  a1þ l Fðx Þ , it follows from (2.2) that

2 2 k k ak bk ak bk k  k k kxk ðak Þ  x k2 6 x  1 þ l Fð~x Þ  x  x  1 þ l Fð~x Þ  x ðak Þ :

ð2:6Þ

And

kxkþ1 ðak Þ  x k2 ¼ ksðxk  x Þ þ ð1  sÞðxk ðak Þ  x Þk2 ¼ s2 kxk  x k2 þ ð1  sÞ2 kxk ðak Þ  x k2 þ 2sð1  sÞhxk  x ; xk ðak Þ  x i: Using the following identity

2ha þ b; bi ¼ ka þ bk2  kak2 þ kbk2 for a ¼ xk  xk ðak Þ; b ¼ xk ðak Þ  x and (2.6), and using 0 < s < 1, we obtain

kxkþ1 ðak Þ  x k2 ¼ s2 kxk  x k2 þ ð1  sÞ2 kxk ðak Þ  x k2 þ sð1  sÞfkxk  x k2  kxk  xk ðak Þk2 þ kxk ðak Þ  x k2 g ¼ skxk  x k2 þ ð1  sÞkxk ðak Þ  x k2  sð1  sÞkxk  xk ðak Þk2

ak bk ~k ak bk ~k 6 skxk  x k2 þ ð1  sÞ kxk  Fðx Þ  x k2  kxk  Fðx Þ  xk ðak Þk2 1þl 1þl  sð1  sÞkxk  xk ðak Þk2

2ak bk k 2ak bk k 6 kxk  x k2  ð1  sÞ kxk  xk ðak Þk2 þ hx ðak Þ  xk ; Fð~xk Þi þ hx  x ; Fð~xk Þi 1þl 1þl



3. The proposed method In this section, we suggest and analyze the new modiﬁed LQP method for solving nonlinear complementarity problems (1.1) and investigate the strategy of how chose the new step size ak . At the kth iteration, LQP method ﬁnds the exact solution for the following system of equations:

bk FðxÞ þ x  ð1  lÞxk  lX 2k x1 ¼ 0:

ð3:1Þ k

We now present an LQP method-based interior prediction-correction method for the NCP. For given x > 0 and bk > 0, each iteration of the proposed method consist of two steps, the ﬁrst step offers a predictor ~ xk and the second step produces the new iterate xkþ1 ðak Þ. 3.1. Prediction step Find an approximate positive solution ~ xk of (3.1), called predictor, such that

0  bk Fð~xk Þ þ ~xk  ð1  lÞxk  lX 2k ð~xk Þ1 ¼: nk

ð3:2Þ

k

and n which satisﬁes

knk k 6 gkxk  ~xk k;

0 < g < 1:

ð3:3Þ

3.2. Correction step For ak > 0 and 0 < s < 1, the new iterate xkþ1 ðak Þ is deﬁned by

 ak bk ~k xkþ1 ðak Þ ¼ sxk þ ð1  sÞPRnþ xk  Fðx Þ : 1þl

ð3:4Þ

Remark 3.1. (3.3) implies that

jhxk  ~xk ; nk ij 6 gkxk  ~xk k2 ;

0 < g < 1:

ð3:5Þ

Remark 3.2. Note that in the special case nk ¼ bk ðFð~ xk Þ  Fðxk ÞÞ, the solution of

bk Fðxk Þ þ ~xk  ð1  lÞxk  lX 2k ð~xk Þ1 ¼ 0;

ð3:6Þ

can be componentwise obtained by

~xkj

¼

ð1  lÞxkj  bk F j ðxk Þ þ

qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ½ð1  lÞxkj  bk F j ðxk Þ2 þ 4lðxkj Þ2 2

:

ð3:7Þ

A. Bnouhachem et al. / Applied Mathematics and Computation 215 (2009) 695–706

699

Moreover for any xk > 0 we have always ~ xk > 0. Now, we begin to investigate how to choose the step size ak in the correction (3.4). For this purpose, we need the following theorem. h i k  2 kþ1 k bk ~k Theorem 3.1. Let xk ðak Þ :¼ P Rnþ xk  a1þ ðak Þ  x k2 , then we have l Fðx Þ and Hðak Þ :¼ kx  x k  kx

Hðak Þ P ð1  sÞWðak Þ P ð1  sÞUðak Þ

ð3:8Þ

where

Wðak Þ ¼ kxk  xk ðak Þk2 þ

2ak bk k hx ðak Þ  ~xk ; Fð~xk Þi; 1þl 

ð3:9Þ

Uðak Þ ¼ 2ak uk  a2k kdk k2 ;

ð3:10Þ

1 1 uk ¼ kxk  ~xk k2 þ hxk  ~xk ; nk i; 1þl 1þl

ð3:11Þ

1 nk : 1þl

ð3:12Þ

and k

d ¼ ðxk  ~xk Þ þ

Proof. By setting q ¼ bk Fð~ xk Þ  nk in (2.3) and y ¼ xk ðak Þ in (2.4), it follows that

hxk ðak Þ  ~xk ;

 1 1 k 1l ðnk  bk Fð~xk ÞÞi 6 kx  xk ðak Þk2  k~xk  xk ðak Þk2  kxk  ~xk k2 : 1þl 2 2ð1 þ lÞ

ð3:13Þ

Recall

hxk ðak Þ  ~xk ; xk  ~xk i ¼

 1 1 k k~x  xk ðak Þk2  kxk  xk ðak Þk2 þ kxk  ~xk k2 : 2 2

ð3:14Þ

Adding (3.13) and (3.14) we then obtain

 xk ðak Þ  ~xk ; xk  ~xk þ

1 ðnk  bk Fð~xk ÞÞ 1þl

 6

l 1þl

kxk  ~xk k2 ;

which implies

 2ak xk ðak Þ  ~xk ; xk  ~xk þ

 1 2ak l k ðnk  bk Fð~xk ÞÞ  kx  ~xk k2 6 0: 1þl 1þl

ð3:15Þ

Using Lemma 2.2 and the deﬁnition of Hðak Þ, we get

Hðak Þ P ð1  sÞ kxk  xk ðak Þk2 þ

 2ak bk k 2ak bk  k x ðak Þ  xk ; Fð~xk Þ þ hx  x ; Fð~xk Þi : 1þl 1þl

ð3:16Þ

Since ~ xk 2 Rnþ and x is a solution of NCP, using the pseudomonotonicity of F we have

h~xk  x ; Fðx Þi ¼ h~xk ; Fðx Þi P 0 ) h~xk  x ; Fð~xk i P 0 and consequently

hxk  x ; Fð~xk Þi P hxk  ~xk ; Fð~xk Þi:

ð3:17Þ

Applying (3.17) to the last term in the right side of (3.16), we obtain



2a b 2 Hðak Þ P ð1  sÞ xk  xk ðak Þ þ k k xk ðak Þ  ~xk ; Fð~xk Þ 1þl



¼ ð1  sÞWðak Þ:

k

Adding Wðak Þ and (3.15) and using the notation of d in (3.12), we get

2

Wðak Þ P xk ðak Þ  xk þ 2ak hxk ðak Þ  xk ; dk i þ 2ak hxk  ~xk ; dk i 

2ak l k kx  ~xk k2 1þl

2ak l k k k k ¼ kxk ðak Þ  xk þ ak d k2  a2k kd k2 þ 2ak hxk  ~xk ; d i  kx  ~xk k2 1þl   1 2ak l k k P 2ak xk  ~xk ; xk  ~xk þ nk  kx  ~xk k2  a2k kd k2 1þl 1þl   1 1 k k ¼ 2ak kxk  ~xk k2 þ hxk  ~xk ; nk i  a2k kd k2 ¼ 2ak uk  a2k kd k2 ¼ Uðak Þ: 1þl 1þl And the conclusion of this theorem is proved. h

ð3:18Þ

700

A. Bnouhachem et al. / Applied Mathematics and Computation 215 (2009) 695–706

Proposition 3.1. Assume that F is continuously differentiable, then we have 2bk k ~k ~k (i) W0 ðak Þ ¼ 1þ l hx ðak Þ  x ; Fðx Þi: (ii) W0 ðak Þ is decreasing function with respect to ak > 0, i.e., Wðak Þ is concave.

xk , let Proof. For given bk ; xk ; ~

 2 a2k b2k ak bk 2ak bk k k k hðak ; yÞ ¼ kFð~xk Þk2  h~x  xk ; Fð~xk Þi y  x  1 þ l Fð~x Þ  2 1þl ð1 þ lÞ

ð3:19Þ

It easy to see that the solution of the following problem

miny fhðak ; yÞjy 2 Rnþ g h i  k bk ~k is y ¼ PRnþ xk  a1þ l Fðx Þ . Substituting y into (3.19) and simplifying it, we have

Wðak Þ ¼ hðak ; yÞj

y¼P Rn

þ

h

i ¼ miny hðak ; yÞjy 2 Rn  þ

a b

k k Fð~ k xk  1þ l x Þ

It follows from  (Chapter 4, Theorem 1.7) that Wðak Þ is differentiable and its derivative is given by

W0 ðak Þ ¼

 @hðak ; yÞ h i @ ak y¼P n xk ak bk Fð~xk Þ R þ

1þl

  2bk ak bk ~k ~k 2ak b2k 2bk xk ðak Þ  xk þ Fðx Þ; Fðx Þ  kFð~xk Þk2  h~xk  xk ; Fð~xk Þi ¼ 1þl 1þl 1þl ð1 þ lÞ2 ¼

2bk hxk ðak Þ  ~xk ; Fð~xk Þi: 1þl 

and the ﬁrst conclusion is proved. We now establish the proof of the second assertion. Let ak > ak > 0, we will prove that

W0 ðak Þ 6 W0 ðak Þ i.e.,

 k  x ðak Þ  xk ðak Þ; Fð~xk Þ 6 0: Setting z :¼ x

k

k bk ~k  a1þ l Fðx Þ;

v :¼

ð3:20Þ

xk ð k Þ

a and z :¼ x

k

k bk ~k  a1þ l Fðx Þ;

  ak bk ~k xk  Fðx Þ  xk ðak Þ; xk ðak Þ  xk ðak Þ 6 0 1þl

v :¼

xk ð k Þ

a in (2.1), respectively, we have ð3:21Þ

and

  ak bk ~k xk  Fðx Þ  xk ðak Þ; xk ðak Þ  xk ðak Þ 6 0: 1þl

ð3:22Þ

Adding (3.21) and (3.22), we obtain

  ðak  ak Þbk xk ðak Þ  xk ðak Þ; xk ðak Þ  xk ðak Þ þ Fð~xk Þ 6 0; 1þl i.e.,

k    x ðak Þ  xk ðak Þ 2 þ ðak  ak Þbk xk ðak Þ  xk ðak Þ; Fð~xk Þ 6 0:     1þl It follows that

 k  x ðak Þ  xk ðak Þ; Fð~xk Þ 6 

1þl xk ðak Þ  xk ðak Þ 2 6 0:  ðak  ak Þbk 

Then, we obtain the inequality (3.20) and complete the proof. h Now for the same kth approximate solution xk , let

ak2 ¼ arg max fUðaÞja > 0g a

ð3:23Þ

A. Bnouhachem et al. / Applied Mathematics and Computation 215 (2009) 695–706

701

and

ak1 ¼ arg maxfWðaÞ j a P 0g:

ð3:24Þ

a

In order to make ak1 be obtained more easily, we approximately compute ak1 by solving the following simple optimization problem:

n

o

ak1 ¼ arg max WðaÞj0 < a 6 m1 ak2 ; a

ð3:25Þ

where m1 P 1. Note that UðaÞ is a quadratic function of a and it reaches its maximum at

ak2 ¼

uk k

kd k2

which is used in½5; 15; 28

ð3:26Þ

and

Uðak2 Þ ¼ ak2 uk :

ð3:27Þ

Based on the Theorem 3.1 and Proposition 3.1, the following convergence results can be proved easily. Proposition 3.2. Let ak1 and ak2 be deﬁned by (3.25) and (3.23), respectively, F be pseudomonotone and continuously differentiable, then we have (i) kxk  x k2  kxkþ1 ðak1 Þ  x k2 P ð1  sÞWðak1 Þ (ii) kxk  x k2  kxkþ1 ðak2 Þ  x k2 P ð1  sÞUðak2 Þ (iii) Wðak1 Þ P Uðak2 Þ Furthermore, if W0 ðak1 Þ ¼ 0, we have (iv) kxk  x k2  kxkþ1 ðak1 Þ  x k2 P ð1  sÞkxk  xkþ1 ðak1 Þk2 . Remark 3.3. The main amount of computation of approximately ﬁnding the point ak1 (at which the maximum value of WðaÞ xk . Note that the value of F at ~ xk will be needed on ð0; m1 ak2  is attained) is to obtain or observe the value of function F at ~ k kþ1 k k   x Þ  Fðx ÞÞ, if we use ak1 instead of ak2 (used in [5,15,21,28]), we don’t again in x ðaÞ. Then, in the special case n ¼ bk ðFð~ need to observe or compute additional value of function F, just cost a relatively minor amount of computation for obtaining projections of some vectors on Rnþ during computing the approximate value of ak1 . According to Proposition 3.2, ak1 can be taken as the more proper step size instead of ak2 . For the convergence analysis of the proposed method, we need the following results. Theorem 3.2 ([5,21]). For given xk 2 Rnþ and bk > 0, let ~ xk and nk satisﬁed to the condition (3.3), then we have the following

ak2 P

1g 2ð1 þ lÞ

ð3:28Þ

and

Uðak2 Þ P

ð1  gÞ2 2ð1 þ lÞ2

kxk  ~xk k2 :

Proof. If hxk  ~ xk ; nk i 6 0, since

l > 0 it follows from (3.3) and (3.12) that

1

k

kd k2 6 kxk  ~xk k2 þ

ð1 þ lÞ2

knk k2 6 kxk  ~xk k2 þ knk k2 6 2kxk  ~xk k2 ;

from (3.26) and (3.30), we obtain

ak2 ¼

uk k

kd k2

P

1g : 2ð1 þ lÞ

Otherwise, if hxk  ~ xk ; nk i P 0, it follows from 0 < l < 1 and (3.3) that

  1 1 1 1 k 1 1 kxk  ~xk k2 þ hxk  ~xk ; nk i P hxk  ~xk ; nk i þ kxk  ~xk k2 kx  ~xk k2 þ 1þl 1þl 1þl 2 1þl 2 ! 1 1 k 1 1 1 k P knk k2 ¼ kx  ~xk k2 þ hxk  ~xk ; nk i þ kd k2 1þl 2 ð1 þ lÞ 2ð1 þ lÞ 2ð1 þ lÞ2

uk ¼

ð3:29Þ

ð3:30Þ

702

A. Bnouhachem et al. / Applied Mathematics and Computation 215 (2009) 695–706

and thus

ak2 P

1 1g P : 2ð1 þ lÞ 2ð1 þ lÞ

It follows from (3.5) and (3.11) that

uk ¼

  1 1 1g kxk  ~xk k2 þ ðxk  ~xk ÞT nk P kxk  ~xk k2 : 1þl 1þl 1þl

Using 3.31, 3.27 and 3.28 directly we obtained (3.29).

ð3:31Þ

h

For fast convergence, we take a relaxation factor c 2 ½1; 2Þ and set the step size ak2 ¼ cak2 . Through simple manipulations we obtain







U cak2 ¼ 2cak2 uk  c2 ak2











ak2 kdk k2 ¼ 2cak2  c2 ak2 uk ¼ cð2  cÞU ak2



ð3:32Þ

It follows from Theorem 3.1, Proposition 3.2, Theorem 3.2 and (3.32) that there is a constant

c :¼

cð2  cÞð1  sÞð1  gÞ2 >0 2ð1 þ lÞ2

such that

kxkþ1  x k2 6 kxk  x k2  ckxk  ~xk k2

8x 2 X :

ð3:33Þ

The following result can be proved by similar arguments as in [5,15,21,28]. Hence the proof is omitted. 1

Theorem 3.3 ([5,15,21,28]). If inf k¼0 bk ¼ b > 0, then the sequence fxk g generated by the proposed method converges to some x1 which is a solution of NCP. The Proposition 3.2 motivates us that we can improve the LQP method by choosing a more proper step size a based on ﬁnding ak1 instead of ak2 which is used in [5,28] and extending the step size by solving the following subproblem:

n

o

ak ¼ max ak1 6 a 6 m2 ak1 j WðaÞ P rWðak1 Þ ;

ð3:34Þ

a

where r 2 ð0; 1Þ and m2 P 2. We now describe the new algorithm as follows. Algorithm 3.1 Step 0. Let b0 > 0, e > 0, 0 < l < 1, 0 < r < 1, 0 < g < 1, 0 < s < 1, m1 P 1, m2 P 2, x0 2 Rnþ and set k :¼ 0. Step 1. If k minðx; FðxÞÞk1 6 , then stop. Otherwise, go to Step 2.  qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ  ~ 2, Step 2. s :¼ ð1  lÞxk  bk Fðxk Þ, xki :¼ si þ ðsi Þ2 þ 4lðxki Þ2 nk :¼ bk ðFð~ xk Þ  Fðxk ÞÞ; while (r > g) bk :¼ bk  0:8=r, s :¼ ð1  lÞxk  bk Fðxk Þ,

r :¼ knk k=kxk  ~ xk k.

~ xki :¼

nk :¼ bk ðFð~ xk Þ  Fðxk ÞÞ;

 qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ  2, si þ ðsi Þ2 þ 4lðxki Þ2

r :¼ knk k=kxk  ~ xk k.

end while Step 3. Searching step size ak :  k ¼ arg maxa fUðaÞ j a > 0g, where UðaÞ is deﬁned by (3.10). Let a Solve the following optimization problem ak ¼ arg maxa fWðaÞ j 0 < a 6 m1 a k g, where WðaÞ is deﬁned by (3.9). Step 4. Extending the step size: ak ¼ maxa fak 6 a 6 m2 ak j WðaÞ P rWðak Þg and

Step 5. bkþ1

k bk ~k xkþ1 ¼ sxk þ ð1  sÞPRnþ ½xk  a1þ l Fðx Þ,

b 0:7 k ; if r 6 0:3; . r ¼ bk ; otherwise:

Step 6. k :¼ k þ 1; go to Step 1.

703

A. Bnouhachem et al. / Applied Mathematics and Computation 215 (2009) 695–706

4. Computational results In this section, we consider two examples to illustrate the efﬁciency and the performance of the proposed algorithm. 4.1. Numerical experiments I We consider the nonlinear complementarity problems:

x P 0;

FðxÞ P 0;

xT FðxÞ ¼ 0;

ð4:1Þ

where

FðxÞ ¼ DðxÞ þ Mx þ q; DðxÞ and Mx þ q are the nonlinear part and linear part of FðxÞ, respectively. The matrix M ¼ AT A þ B is computed as follows. A is n  n matrix whose entries are randomly generated in the interval (5, +5) and the skew-symmetric matrix B is generated in the same way. The components of D(x) are Dj ðxÞ ¼ dj  arctanðxj Þ and dj is chosen randomly in (0,1). In all tests we take the logarithmic proximal parameter l ¼ 0:1; s ¼ 0:01 and all iterations start with x0 ¼ ð1; . . . ; 1ÞT and b0 ¼ 1. The stopping criterion was set to be

k minfx; FðxÞgk1 6 107 : k minfx0 ; Fðx0 Þgk1 All codes are written in Matlab and run on a P4-2.00G note book computer. We test the problems (4.1) with different dimensions and q 2 ð500; 500Þ in Table 1, and q 2 ð500; 0Þ in Table 2. We compared the proposed method with that of Xu et al. method . In all tests we take r ¼ 0:05; m1 ¼ 3; m2 ¼ 4; g ¼ 0:9, c ¼ 1:8. The test results for problems (4.1) are reported in Tables 1 and 2. k is the number of iterations and l denotes the number of evaluations of mapping F. 4.2. Numerical experiments II In this subsection, we apply the proposed method to the trafﬁc equilibrium problems and present corresponding numerical results. Consider a network [N; L] of nodes N and directed links L, which consists of a ﬁnite sequence of connecting links with a certain orientation. Let a; b, etc., denote the links, and let p; q, etc., denote the paths. We let x denote an origin/destination (O/D) pair of nodes of the network and Px denotes the set of all paths connecting O/D pair x. Note that the path-arc incidence matrix and the path-O/D pair incidence matrix, denoted by A and B, respectively, are determined by the given network and O/D pairs. To see how to convert a trafﬁc equilibrium problem into a variational inequality, we take into account a simple example depicted in Fig. 1.

Table 1 The numerical results for problem (4.1) with q 2 ð500; 500Þ. n

200 300 400 500 700 1000

The proposed method

The method in 

k

l

k

l

223 247 248 272 261 249

486 538 540 594 570 544

385 416 418 468 442 423

818 884 890 997 941 899

Table 2 The numerical results for problem (4.1) with q 2 ð500; 0Þ. n

200 300 400 500 700 1000

The proposed method

The method in 

k

l

k

l

428 422 518 552 483 514

930 916 1124 1200 1051 1118

724 701 857 923 795 856

1541 1492 1823 1962 1690 1818

704

A. Bnouhachem et al. / Applied Mathematics and Computation 215 (2009) 695–706

Fig. 1. An illustrative example of given directed network and the O/D pairs.

For the given example in Fig. 1, the path-arc incidence matrix A and the path-O/D pair incidence matrix B have the following forms:

Let xp represent the trafﬁc ﬂow on path p and fa denote the link load on link a, then the arc-ﬂow vector f is given by

f ¼ AT x:

ð4:2Þ

Let dx denote the trafﬁc amount between O/D pair x, which must satisfy

dx ¼

X

xp :

ð4:3Þ

p2P x

Thus, the O/D pair-trafﬁc amount vector d is given by

d ¼ BT x:

ð4:4Þ

Let tðf Þ ¼ ft a ; a 2 Lg be the vector of link travel costs, which is a function of the link ﬂow. A user traveling on path p incurs a (path) travel cost hp . For given link travel cost vector t, the path travel cost vector h is given by

h ¼ Atðf Þ and thus hðxÞ ¼ AtðAT xÞ:

ð4:5Þ

Associated with every O/D pair x, there is a travel disutility kx ðdÞ. Since both the path costs and the travel disutilities are functions of the ﬂow pattern x, the trafﬁc network equilibrium problem is to seek the path ﬂow pattern x such that

x P 0 ðx  x ÞT Fðx Þ P 0;

8x P 0

ð4:6Þ

where

F p ðxÞ ¼ hp ðxÞ  kx ðdðxÞÞ;

8 x;

p 2 Px

and thus

Fig. 2. A directed network with 25 nodes and 37 links.

ð4:7Þ

705

A. Bnouhachem et al. / Applied Mathematics and Computation 215 (2009) 695–706

FðxÞ ¼ AtðAT xÞ  BkðBT xÞ:

ð4:8Þ

We apply the proposed method to the example taken from  (Example 7.5 in ), which consisted of 25 nodes, 37 links and 6 O/D pairs. The network is depicted in Fig. 2. For this example, there are together 55 paths for the 6 given O/D pairs and hence the dimension of the variable x is 55. Therefore, the path-arc incidence matrix A is a 55  37 matrix and the path-O/D pair incidence matrix B is a 55  6 matrix. The user cost of traversing link a is given in Table 3. The disutility function is given by

kx ðdÞ ¼ mx dx þ qx

ð4:9Þ

and the coefﬁcients mx and qx in the disutility function of different O/D pairs for this example are given in Table 4. The test results for problems (4.6) for different e are reported in Table 5, k is the number of iterations and l denotes the number of evaluations of mapping F. The stopping criterion is

k minfx; FðxÞgk1 6 e: k minfx0 ; Fðx0 Þgk1 Tables 1and 2 show that the proposed method is very efﬁcient algorithm even for large-scale classical NCP. Moreover, the new step size ak plays important role to reduce the iterative numbers and the evaluation numbers of F. Table 5 shows that the new method is more ﬂexible and efﬁcient to solve trafﬁc equilibrium problem. Moreover, it demonstrates computationally that the new method is more effective than the method presented in  in the sense that the new method needs fewer iteration and less evaluation numbers of F, which clearly illustrate its efﬁciency.

Table 3 The link traversing cost functions ta ðf Þ in the example. t1 ðf Þ ¼ 5  105 f14 þ 5f 1 þ 2f 2 þ 500 t2 ðf Þ ¼ 3  105 f24 þ 4f 2 þ 4f 1 þ 200 t3 ðf Þ ¼ 5  105 f34 þ 3f 3 þ f4 þ 350 t4 ðf Þ ¼ 3  105 f44 þ 6f 4 þ 3f 5 þ 400 t5 ðf Þ ¼ 6  105 f54 þ 6f 5 þ 4f 6 þ 600 t6 ðf Þ ¼ 7f 6 þ 3f 7 þ 500 t7 ðf Þ ¼ 8  105 f74 þ 8f 7 þ 2f 8 þ 400 t8 ðf Þ ¼ 4  105 f84 þ 5f 8 þ 2f 9 þ 650 t9 ðf Þ ¼ 105 f94 þ 6f 9 þ 2f 1 0 þ 700 t10 ðf Þ ¼ 4f 10 þ f12 þ 800 4 þ 7f 11 þ 4f 12 þ 650 t11 ðf Þ ¼ 7  105 f11 t12 ðf Þ ¼ 8f 12 þ 2f 13 þ 700 5 4 t13 ðf Þ ¼ 10 f13 þ 7f 13 þ 3f 18 þ 600 t14 ðf Þ ¼ 8f 14 þ 3f 15 þ 500 4 þ 9f 15 þ 2f 14 þ 200 t15 ðf Þ ¼ 3  105 f15 t16 ðf Þ ¼ 8f 16 þ 5f 12 þ 300 4 þ 7f 17 þ 2f 15 þ 450 t17 ðf Þ ¼ 3  105 f17 t18 ðf Þ ¼ 5f 18 þ f16 þ 300 t19 ðf Þ ¼ 8f 19 þ 3f 17 þ 600

t 20 ðf Þ ¼ 3 t 21 ðf Þ ¼ 4 t 22 ðf Þ ¼ 2 t 23 ðf Þ ¼ 3 t 24 ðf Þ ¼ 2 t 25 ðf Þ ¼ 3 t 26 ðf Þ ¼ 6 t 27 ðf Þ ¼ 3 t 28 ðf Þ ¼ 3 t 29 ðf Þ ¼ 3 t 30 ðf Þ ¼ 4 t 31 ðf Þ ¼ 3 t 32 ðf Þ ¼ 6 t 33 ðf Þ ¼ 4 t 34 ðf Þ ¼ 6 t 35 ðf Þ ¼ 3 t 36 ðf Þ ¼ 2 t 37 ðf Þ ¼ 6

                 

4 105 f20 þ 6f 20 þ f21 þ 300 4 105 f21 þ 4f 21 þ f22 þ 400 4 105 f22 þ 6f 22 þ f23 þ 500 5 4 10 f23 þ 9f 23 þ 2f 24 þ 350 5 4 10 f24 þ 8f 24 þ f25 þ 400 4 105 f25 þ 9f 25 þ 3f 26 þ 450 4 105 f26 þ 7f 26 þ 8f 27 þ 300 4 105 f27 þ 8f 27 þ 3f 28 þ 500 4 105 f28 þ 7f 28 þ 650 4 105 f29 þ 3f 29 þ f30 þ 450 5 4 10 f30 þ 7f 30 þ 2f 31 þ 600 5 4 10 f31 þ 8f 31 þ f32 þ 750 4 105 f32 þ 8f 32 þ 3f 33 þ 650 4 105 f33 þ 9f 33 þ 2f 31 þ 750 4 105 f34 þ 7f 34 þ 3f 30 þ 550 4 105 f35 þ 8f 35 þ 3f 32 þ 600 4 105 f36 þ 8f 36 þ 4f 31 þ 750 5 4 10 f37 þ 5f 37 þ f36 þ 350

Table 4 The O/D pairs and the parameters in (4.9) of the example. (O,D) Pair x

(1,20)

(1,25)

(2,20)

(3,25)

(1,24)

(11,25)

mx qx jP x j

1 1000 10

6 800 15

10 2000 9

5 6000 6

7 8000 10

9 7000 5

Table 5 Numerical results for different e. Different

5

10 106 107 108 109

e

The proposed method

The method in 

k

l

k

l

236 306 376 446 516

523 677 831 985 1139

401 519 643 757 877

864 1119 1386 1632 1892

706

A. Bnouhachem et al. / Applied Mathematics and Computation 215 (2009) 695–706

5. Conclusions In this paper, we propose a new modiﬁed logarithmic-quadratic proximal (LQP) method for solving nonlinear complementarity problems (NCP) by using a new step size. In order to obtain the new step size, we do not need to compute additionally the value of function F. This is very important especially in practical problems in which the cost of computing or observing the value of function F is very expensive. Furthermore, numerical experiments show that the proposed method is attractive in practice. How to design other efﬁcient LQP methods for solving NCP is worthy of further investigations in the future. Acknowledgements This research was supported partly by NSFC Grants Nos: 70731002 and 70571034, ‘‘The Second Term’985’ EngineeringInnovation Center for Computational Laboratory of Evolutionary Economic Theory and Its Applications” and ‘‘Study on the Evolution of Complex Economic System” at ‘‘Innovation Center of Economic Transition and Development of Nanjing University” of Ministry of Education, China. The authors would like to express their gratitude to the referees for their valuable suggestions and constructive comments. References  A. Auslender, Optimization méthodes numériques, Mason, Paris, 1976.  A. Auslender, M. Teboulle, S. Ben-Tiba, A logarithmic-quadratic proximal method for variational inequalities, Comput. Optim. Appl. 12 (1999) 31–40.  A. Auslender, M. Teboulle, S. Ben-Tiba, Interior proximal and multiplier methods based on second order homogenous kernels, Math. Oper. Res. 24 (1999) 646–668.  A. Auslender, M. Teboule, Interior projection-like methods for monotone variational inequalities, Math. Prog. Ser. A 104 (2005) 39–68.  A. Bnouhachem, An LQP Method for pseudomonotone variational inequalities, J. Global Optim. 36 (3) (2006) 351–363.  A. Bnouhachem, X. Yuan, An extended LQP method for monotone nonlinear complementarity problems, J. Optim. Theory Appl. 135 (3) (2007) 343–353.  R.S. Burachik, A.N. Iusem, B.F. Svaiter, Enlargement of monotone operators with applications to variational inequalities, Set-Valued Anal. 5 (1997) 159– 180.  R.S. Burachik, A.N. Iusem, A generalized proximal point algorithm for the variational inequality problem in a Hilbert space, SIAM J. Optim. 8 (1998) 197–216.  J. Eckestein, Approximate iterations in Bregman function-based proximal algorithms, Math. Prog. 83 (1998) 113–123.  M.C. Ferris, J.S. Pang, Engineering and economic applications of complementary problems, SIAM Rev. 39 (1997) 669–713.  A. Fischer, Solution of monotone complementarity problems with locally Lipschitzian functions, Math. Prog. 76 (1997) 513–532.  R. Glowinski, Numerical methods for nonlinear variational inequality problems, Springer-Verlag, New York, NY, 1984.  O. Guler, On the convergence of the proximal point algorithm for convex minimization, SIAM J. Control. Optim. 29 (1991) 403–419.  D.R. Han, B.S. He, A new accuracy criterion for approximate proximal point algorithms, J. Math. Anal. Appl. 263 (2001) 343–354.  B.S. He, L.-Z. Liao, X.M. Yuan, A LQP based interior prediction-correction method for nonlinear complementarity problems, J. Comput. Math. 24 (1) (2006) 33–44.  P.T. Harker, J.S. Pang, Finite-dimensional variational inequality and nonlinear complementarity problems: a survey of theory, algorithms and applications, Math. Prog. 48 (1990) 161–220.  J.J. Moré, Global methods for nonlinear complementarity problems, Math. Oper. Res. 21 (1996) 589–614.  B. Martinet, determination approachée d’un point ﬁxe d’une application pseudo-contractante, C.R. Acad. Sci. Paris 274 (1972) 163–165.  M.A. Noor, General variational inequalities, Appl. Math. Lett. 1 (1988) 119–121.  M.A. Noor, Some developments in general variational inequalities, Appl. Math. Comput. 152 (2004) 199–277.  M.A. Noor, A. Bnouhachem, Modiﬁed proximal point method for nonlinear complementarity problems, J. Comput. Appl. Math. 197 (2) (2006) 395–405.  M.A. Noor, K.I. Noor, T.M. Rassias, Some aspects of variational inequalities, J. Comput. Appl. Math. 47 (1993) 285–312.  R.T. Rockafellar, Monotone operators and the proximal point algorithm, SIAM J. Control. Optim. 14 (1976) 877–898.  P.J. Da Silva e Silva, J. Eckstein, C. Humes, Rescaling and stepsize selection in proximal methods using separable generalized distance, Rutcor Research Report Rrr 35-99, Rutgers University, 1999.  M.V. Solodov, B.F. Svaiter, An inexact hybrid general proximal point algorithm and some new result on the theory of Bregman functions, Math. Oper. Res. 25 (2000) 214–230.  M.V. Solodov, B.F. Svaiter, Error bounds for proximal point subproblems and associated inexact proximal point algorithms, Math. Prog. Ser. B 88 (2000) 371–389.  M. Teboulle, Convergence of proximal-like algorithms, SIAM J. Optim. 7 (1997) 1069–1083.  Y. Xu, B.S. He, X. Yuan, A hybrid inexact logarithmic-quadratic proximal method for nonlinear complementarity problems, J. Math. Anal. Appl. 322 (2006) 276–287.  M.H. Xu, X.M. Yuan, Q.L. Huang, An improved general extra-gradient method with reﬁned step size for nonlinear monotone variational inequalities, J. Global Optim. 39 (2007) 155–169.