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 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 modified 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 significantly 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 efficiency of the proposed method. This new method can be considered as a significant refinement 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 finding the zero of an appropriately maximal monotone operator T ¼ F þ N Rnþ , namely, find x 2 Rnþ such that 0 2 Tðx Þ, where N Rnþ ð:Þ is the normal cone operator to Rnþ defined by

NRnþ ðxÞ :¼

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



8v 2 Rnþ ; if

x 2 Rnþ ;

otherwise:

A well known method to find the zero of a maximal monotone operator T is the proximal point algorithm (PPA). It was first introduced by Martinet [18] and further extended by Rockafellar [23] 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 [9] 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) [23] 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 [23] proposed the first relative error criterion of PPA that all the relaxation factors are summable. Eckstein [9] 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 [14] proved the convergence under the following accuracy criterion

kek k 6 gk kxk  xk1 k with

þ1 X

g2k < þ1:

ð1:6Þ

k¼0

Furthermore, Da Silva et al. [24] 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. [3] have proposed a new type of proximal interior method through replacing the quadratic function (1.3) by d/ ðx; xk Þ which could be defined as

d/ ðx; yÞ ¼

n X

y2j /ðy1 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 efficient algorithms to solve NCP. Let m > l > 0 be given fixed parameters, and define

(

m ðt

2

/ðtÞ ¼

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

þ1

t>0

otherwise:

In [2], 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 defined 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 x1 ;  where X k ¼ diag xk1 ; xk2 ; . . . ; xkn , x1 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 x1 ¼ 0:

ð1:9Þ

In Auslender, Teboulle and Ben-Tiba’s original work [2], it is required that the subproblem should be solved exactly. To release the difficulty in Auslender’s problem, very recently, He et al. [15], Bnouhachem [5], Noor and Bnouhachem [21],

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

697

and Xu et al. [28] introduced some LQP based prediction-correction methods which do not suffer from the above difficulty 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 significantly relaxed accuracy criterion and the new iterate is computed directly by an explicit formula derived from the original LQP method for [15], 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 significant refinement and improvement of the method of [28]. We also study the global convergence of the proposed modified LQP method under some mild conditions. Some preliminary computational results are given to illustrate the efficiency 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 first 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Þ

Definition 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 [2]. 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 x1 ¼ 0;

ð2:3Þ

then for any y P 0 we have

hy  x; qi P

 1l 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 defined 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 modified 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 finds the exact solution for the following system of equations:

bk FðxÞ þ x  ð1  lÞxk  lX 2k x1 ¼ 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 first 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 satisfies

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 defined 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 Þ þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ½ð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 1l ð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 definition 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 [1] (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 first conclusion is proved. We now establish the proof of the second assertion. Let ak > ak > 0, we will prove that

W0 ðak Þ 6 W0 ðak Þ i.e.,

 k  x ðak Þ  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 Þ;

  ak bk ~k xk  Fðx Þ  xk ðak Þ; xk ðak Þ  xk ðak Þ 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 ðak Þ  xk ðak Þ 6 0: 1þl

ð3:22Þ

Adding (3.21) and (3.22), we obtain

  ðak  ak Þbk xk ðak Þ  xk ðak Þ; xk ðak Þ  xk ðak Þ þ Fð~xk Þ 6 0; 1þl i.e.,

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

 k  x ðak Þ  xk ðak Þ; Fð~xk Þ 6 

1þl xk ðak Þ  xk ðak Þ 2 6 0:  ðak  ak Þbk 

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

ak2 ¼ arg max fUðaÞja > 0g a

ð3:23Þ

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

701

and

ak1 ¼ arg maxfWðaÞ j a P 0g:

ð3:24Þ

a

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

n

o

ak1 ¼ arg max WðaÞj0 < a 6 m1 ak2 ; a

ð3:25Þ

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

ak2 ¼

uk k

kd k2

which is used in½5; 15; 28

ð3:26Þ

and

Uðak2 Þ ¼ ak2 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 ak1 and ak2 be defined by (3.25) and (3.23), respectively, F be pseudomonotone and continuously differentiable, then we have (i) kxk  x k2  kxkþ1 ðak1 Þ  x k2 P ð1  sÞWðak1 Þ (ii) kxk  x k2  kxkþ1 ðak2 Þ  x k2 P ð1  sÞUðak2 Þ (iii) Wðak1 Þ P Uðak2 Þ Furthermore, if W0 ðak1 Þ ¼ 0, we have (iv) kxk  x k2  kxkþ1 ðak1 Þ  x k2 P ð1  sÞkxk  xkþ1 ðak1 Þk2 . Remark 3.3. The main amount of computation of approximately finding the point ak1 (at which the maximum value of WðaÞ xk . Note that the value of F at ~ xk will be needed on ð0; m1 ak2  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 ak1 . According to Proposition 3.2, ak1 can be taken as the more proper step size instead of ak2 . 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 satisfied to the condition (3.3), then we have the following

ak2 P

1g 2ð1 þ lÞ

ð3:28Þ

and

Uðak2 Þ 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

ak2 ¼

uk k

kd k2

P

1g : 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

ak2 P

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

It follows from (3.5) and (3.11) that

uk ¼

  1 1 1g 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 ak2 ¼ cak2 . Through simple manipulations we obtain







U cak2 ¼ 2cak2 uk  c2 ak2











ak2 kdk k2 ¼ 2cak2  c2 ak2 uk ¼ cð2  cÞU ak2



ð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 finding ak1 instead of ak2 which is used in [5,28] and extending the step size by solving the following subproblem:

n

o

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

ð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.  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ~ 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 ÞÞ;

 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2, si þ ðsi Þ2 þ 4lðxki Þ2

r :¼ knk k=kxk  ~ xk k.

end while Step 3. Searching step size ak :  k ¼ arg maxa fUðaÞ j a > 0g, where UðaÞ is defined by (3.10). Let a Solve the following optimization problem ak ¼ arg maxa fWðaÞ j 0 < a 6 m1 a k g, where WðaÞ is defined by (3.9). Step 4. Extending the step size: ak ¼ maxa fak 6 a 6 m2 ak j WðaÞ P rWðak Þ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 efficiency 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 107 : 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 [28]. 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 traffic equilibrium problems and present corresponding numerical results. Consider a network [N; L] of nodes N and directed links L, which consists of a finite 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 traffic 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 [28]

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 [28]

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 traffic flow on path p and fa denote the link load on link a, then the arc-flow vector f is given by

f ¼ AT x:

ð4:2Þ

Let dx denote the traffic amount between O/D pair x, which must satisfy

dx ¼

X

xp :

ð4:3Þ

p2P x

Thus, the O/D pair-traffic 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 flow. 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 flow pattern x, the traffic network equilibrium problem is to seek the path flow 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 [19] (Example 7.5 in [19]), 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 coefficients 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 efficient 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 flexible and efficient to solve traffic equilibrium problem. Moreover, it demonstrates computationally that the new method is more effective than the method presented in [28] in the sense that the new method needs fewer iteration and less evaluation numbers of F, which clearly illustrate its efficiency.

Table 3 The link traversing cost functions ta ðf Þ in the example. t1 ðf Þ ¼ 5  105 f14 þ 5f 1 þ 2f 2 þ 500 t2 ðf Þ ¼ 3  105 f24 þ 4f 2 þ 4f 1 þ 200 t3 ðf Þ ¼ 5  105 f34 þ 3f 3 þ f4 þ 350 t4 ðf Þ ¼ 3  105 f44 þ 6f 4 þ 3f 5 þ 400 t5 ðf Þ ¼ 6  105 f54 þ 6f 5 þ 4f 6 þ 600 t6 ðf Þ ¼ 7f 6 þ 3f 7 þ 500 t7 ðf Þ ¼ 8  105 f74 þ 8f 7 þ 2f 8 þ 400 t8 ðf Þ ¼ 4  105 f84 þ 5f 8 þ 2f 9 þ 650 t9 ðf Þ ¼ 105 f94 þ 6f 9 þ 2f 1 0 þ 700 t10 ðf Þ ¼ 4f 10 þ f12 þ 800 4 þ 7f 11 þ 4f 12 þ 650 t11 ðf Þ ¼ 7  105 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  105 f15 t16 ðf Þ ¼ 8f 16 þ 5f 12 þ 300 4 þ 7f 17 þ 2f 15 þ 450 t17 ðf Þ ¼ 3  105 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 105 f20 þ 6f 20 þ f21 þ 300 4 105 f21 þ 4f 21 þ f22 þ 400 4 105 f22 þ 6f 22 þ f23 þ 500 5 4 10 f23 þ 9f 23 þ 2f 24 þ 350 5 4 10 f24 þ 8f 24 þ f25 þ 400 4 105 f25 þ 9f 25 þ 3f 26 þ 450 4 105 f26 þ 7f 26 þ 8f 27 þ 300 4 105 f27 þ 8f 27 þ 3f 28 þ 500 4 105 f28 þ 7f 28 þ 650 4 105 f29 þ 3f 29 þ f30 þ 450 5 4 10 f30 þ 7f 30 þ 2f 31 þ 600 5 4 10 f31 þ 8f 31 þ f32 þ 750 4 105 f32 þ 8f 32 þ 3f 33 þ 650 4 105 f33 þ 9f 33 þ 2f 31 þ 750 4 105 f34 þ 7f 34 þ 3f 30 þ 550 4 105 f35 þ 8f 35 þ 3f 32 þ 600 4 105 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 106 107 108 109

e

The proposed method

The method in [28]

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 modified 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 efficient 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 [1] A. Auslender, Optimization méthodes numériques, Mason, Paris, 1976. [2] A. Auslender, M. Teboulle, S. Ben-Tiba, A logarithmic-quadratic proximal method for variational inequalities, Comput. Optim. Appl. 12 (1999) 31–40. [3] 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. [4] A. Auslender, M. Teboule, Interior projection-like methods for monotone variational inequalities, Math. Prog. Ser. A 104 (2005) 39–68. [5] A. Bnouhachem, An LQP Method for pseudomonotone variational inequalities, J. Global Optim. 36 (3) (2006) 351–363. [6] A. Bnouhachem, X. Yuan, An extended LQP method for monotone nonlinear complementarity problems, J. Optim. Theory Appl. 135 (3) (2007) 343–353. [7] 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. [8] 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. [9] J. Eckestein, Approximate iterations in Bregman function-based proximal algorithms, Math. Prog. 83 (1998) 113–123. [10] M.C. Ferris, J.S. Pang, Engineering and economic applications of complementary problems, SIAM Rev. 39 (1997) 669–713. [11] A. Fischer, Solution of monotone complementarity problems with locally Lipschitzian functions, Math. Prog. 76 (1997) 513–532. [12] R. Glowinski, Numerical methods for nonlinear variational inequality problems, Springer-Verlag, New York, NY, 1984. [13] O. Guler, On the convergence of the proximal point algorithm for convex minimization, SIAM J. Control. Optim. 29 (1991) 403–419. [14] D.R. Han, B.S. He, A new accuracy criterion for approximate proximal point algorithms, J. Math. Anal. Appl. 263 (2001) 343–354. [15] 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. [16] 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. [17] J.J. Moré, Global methods for nonlinear complementarity problems, Math. Oper. Res. 21 (1996) 589–614. [18] B. Martinet, determination approachée d’un point fixe d’une application pseudo-contractante, C.R. Acad. Sci. Paris 274 (1972) 163–165. [19] M.A. Noor, General variational inequalities, Appl. Math. Lett. 1 (1988) 119–121. [20] M.A. Noor, Some developments in general variational inequalities, Appl. Math. Comput. 152 (2004) 199–277. [21] M.A. Noor, A. Bnouhachem, Modified proximal point method for nonlinear complementarity problems, J. Comput. Appl. Math. 197 (2) (2006) 395–405. [22] M.A. Noor, K.I. Noor, T.M. Rassias, Some aspects of variational inequalities, J. Comput. Appl. Math. 47 (1993) 285–312. [23] R.T. Rockafellar, Monotone operators and the proximal point algorithm, SIAM J. Control. Optim. 14 (1976) 877–898. [24] 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. [25] 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. [26] 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. [27] M. Teboulle, Convergence of proximal-like algorithms, SIAM J. Optim. 7 (1997) 1069–1083. [28] 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. [29] M.H. Xu, X.M. Yuan, Q.L. Huang, An improved general extra-gradient method with refined step size for nonlinear monotone variational inequalities, J. Global Optim. 39 (2007) 155–169.