www.elsevier.com/locate/cam
A regularization semismooth Newton method based on the generalized Fischer–Burmeister function for P 0 -NCPs
Jein-Shan Chen
a,∗,1, Shaohua Pan
baDepartment of Mathematics, National Taiwan Normal University, Taipei 11677, Taiwan bSchool of Mathematical Sciences, South China University of Technology, Guangzhou 510641, China
Received 1 June 2007; received in revised form 29 August 2007
Abstract
We consider a regularization method for nonlinear complementarity problems with F being a P0-function which replaces the original problem with a sequence of the regularized complementarity problems. In this paper, this sequence of regularized comple- mentarity problems are solved approximately by applying the generalized Newton method for an equivalent augmented system of equations, constructed by the generalized Fischer–Burmeister (FB) NCP-functionspwith p > 1. We test the performance of the regularization semismooth Newton method based on the family of NCP-functions through solving all test problems from MCPLIB.
Numerical experiments indicate that the method associated with a smaller p, for example p ∈ [1.1, 2], usually has better numerical performance, and the generalized FB functionspwith p ∈ [1.1, 2) can be used as the substitutions for the FB function2.
© 2007 Elsevier B.V. All rights reserved.
Keywords: Nonlinear complementarity problem (NCP); Generalized Fischer–Burmeister function; P0-function; Semismooth Newton method
1. Introduction
The nonlinear complementarity problem (NCP) is to find a point x ∈ Rnsuch that
x 0, F (x)0, x, F (x) = 0, (1)
where·, · is the Euclidean inner product and F = (F1, F2, . . . , Fn)Tis a map fromRntoRn. We assume that F is continuously differentiable throughout this paper. The NCP has attracted much attention due to its various applications in operations research, economics, and engineering[12,17,24]. There have been many methods proposed for solving the NCP, including merit function approaches[16,21,23,33], nonsmooth Newton methods[11,22,34], smoothing methods [5,18,27,32]and regularization methods[9,19,29,30]. All the aforementioned methods usually exploit so-called NCP- functions defined as below.
Definition 1.1. A function : R2 → R is called an NCP-function (or C-function standing for Complementarity function, see[10]) if it satisfies
(a, b) = 0 ⇐⇒ a 0, b0, ab = 0. (2)
∗Corresponding author.
E-mail addresses:[email protected](J. Chen),[email protected](S. Pan).
1Member of Mathematics Division, National Center for Theoretical Sciences, Taipei Office.
0377-0427/$ - see front matter © 2007 Elsevier B.V. All rights reserved.
doi:10.1016/j.cam.2007.08.020
Over the past two decades, a variety of NCP-functions has been studied; see[15,31]and references therein. Among which, a popular NCP-function intensively studied is the well-known Fischer–Burmeister (FB) NCP-function[13,14]
defined as
FB(a, b) =
a2+ b2− (a + b). (3)
SinceFBsatisfies (2), the NCP is equivalent to a system of nonsmooth equations
FB(x) :=
⎛
⎜⎜
⎜⎝
FB(x1, F1(x))
··
FB(xn·, Fn(x))
⎞
⎟⎟
⎟⎠= 0. (4)
Then we have the merit functionFB: Rn→ R+for the NCP, defined by
FB(x) := 1
2 FB(x) 2=1 2
n i=1
FB(xi, Fi(x))2. (5)
Recently, a family of new NCP-functions based on the FB function (3) were studied in[2,6]. In particular, they define
p: R2→ R by
p(a, b) := (a, b) p− (a + b), (6)
where p is any fixed real number in the interval (1, +∞) and (a, b) pdenotes the p-norm of (a, b), namely, (a, b) p=
√p
|a|p+ |b|p. In other words, in the functionp, we replace the Euclidean norm of (a, b) in the FB function (3) by a more general p-norm with p ∈ (1, +∞). Similarly, the NCP is equivalent to the nonsmooth system
p(x) :=
⎛
⎜⎜
⎜⎝
p(x1, F1(x))
··
p(xn, F· n(x))
⎞
⎟⎟
⎟⎠= 0, (7)
which induces a family of merit functionsp : Rn→ R for the NCP as below
p(x) := 1
2 p(x) 2=1 2
n i=1
p(xi, Fi(x))2. (8)
As seen in[6], the merit functionpfor any given p > 1 enjoys all favorable properties as the FB merit function FB
holds. Moreover, numerical experiments there indicate that the descent method based on the merit functionp has better performance when p decreases in (1, +∞). However, it is still unknown whether such phenomenon occurs in other approaches for the NCP. The main purpose of this paper is to investigate how the generalized FB NCP-functions
pwith p ∈ (1, +∞) behave in a regularization semismooth Newton method for solving the NCP.
It is well known that the regularization approach is designed to handle ill-posed problems which substitutes the solution of original problem with the solution of a sequence of well-posed problems whose solutions converging to the solution of the original problem; see[4,3,9,19,30]and references therein. In the context of complementarity problems, if we consider the so-called Tikhonov regularization, this scheme consists of solving a sequence of complementarity problems NCP(Fε):
x 0, Fε(x)0, x, Fε(x) = 0, (9)
where ε > 0 is a parameter tending to zero and Fεis given by
Fε(x) := F (x) + εx. (10)
Let Fε,i(x) denote the ith component of Fε(x) and define the map p,ε: Rn→ Rnby
p,ε(x) :=
⎛
⎜⎜
⎜⎝
p(x1, Fε,1(x))
··
p(xn, F·ε,n(x))
⎞
⎟⎟
⎟⎠. (11)
Then the regularized problem NCP(Fε) for any given ε > 0 can be reformulated as
p,ε(x) = 0,
which leads to a merit functionp,ε: Rn→ R+for the NCP(Fε):
p,ε(x) := 1
2 p,ε(x) 2=1 2
n i=1
p(xi, Fε,i(x))2. (12)
Therefore, the original NCP is actually equivalent to solving a sequence of nonsmooth systems of equationsp,ε(x)=0 with ε approaching to 0. From this, we see that the parameter ε plays the same role as the smoothing parameter in smoothing methods for the NCP, except that ε is imposed on the mapping F instead of the NCP-function p.
In this paper, the sequence of subproblemsp,ε(x) = 0 with ε tending to 0 will be solved approximately by applying the generalized Newton method for an augmented system of equations equivalent to the NCP. Specifically, we let z := (ε, x) ∈ R+× Rnby viewing ε as a variable, and define the mapping Hp: R+× Rn→ Rn+1by
Hp(z) :=
⎡
⎢⎢
⎣
p(x1, Fεε,1(x)) ...
p(xn, Fε,n(x))
⎤
⎥⎥
⎦ . (13)
Notice that if the function p,ε(x) defined by (11) is viewed as a function of ε and x, then we may denote it as
p(z) : =p(ε, x) = p,ε(x). Hence, (13) is the same as
Hp(z) =
ε
p(z)
.
It is easily verified that the NCP is equivalent to the augmented system of equations
Hp(z) = Hp(ε, x) = 0, (14)
which naturally induces a merit function Gp: Rn+1→ R+given by Gp(z) =1
2 Hp(z) 2=1
2(ε2+ p,ε(x) 2) =1
2ε2+ p(z). (15)
The function Hp is locally Lipschitz continuous sincep is locally Lipschitz continuous (see[6]). Furthermore, as shown in Section 3, Hpis semismooth. By this, we apply the generalized Newton method developed by[26,28]for (14), and establish a regularized semismooth Newton-type algorithm which in each step solves a regularized problem NCP(Fε) approximately. Compared with the semismooth Newton method based on (7), the method has a remarkable advantage in handling the P0-NCPs (see Section 4) since the merit functionp,ε(x) has bounded level sets for such NCPs. We examine the numerical performance of the algorithm by applying it for all test problems from MCPLIB with three specific NCP-functions1.1, 2 and5. Numerical results indicate that the method associated with a smaller p, for example p ∈ [1.1, 2], usually has better numerical performance, and the generalized FB functions p with p ∈ [1.1, 2) can be used as the substitutions for the FB function 2.
Throughout this paper,R+andR++denote the set of nonnegative real numbers and the set of positive real numbers, respectively; Rn represents the space of n-dimensional real column vectors; and T is the transpose notation. For any differentiable function f : Rn → R, ∇f (x) denotes the gradient of f at x. For any differentiable mapping
F = (F1, . . . , Fm)T : Rn → Rm, F (x) means the Jacobian matrix of F at x while ∇F (x) = [∇F1(x) . . . ∇Fm(x)]
denotes the transpose Jacobian of F at x. If W is an n × n matrix with entries Wj k, j, k = 1, 2, . . . , n, and J and K are index sets such that J, K ⊆ {1, 2, . . . , n}, we denote by WJKthe|J| × |K| submatrix of W consisting of entries Wj k, j ∈ J, k ∈ K. We denote by x p the p-norm of x and by x the Euclidean norm of x. In addition, unless otherwise stated, we always assume that p is any fixed real number in (1, +∞) and denote S∗by the solution set of the NCP if it is nonempty.
2. Preliminaries
In this section, we recall some background concepts and materials which will be used in the subsequent analysis.
We start with the definition of P -matrix and P0-matrix.
Definition 2.1. Given a matrix M ∈ Rn×n, then M is a (a) P0-matrix if each of its principal minors is nonnegative;
(b) P -matrix if each of its principal minors is positive.
Clearly, a positive semidefinite matrix is a P0-matrix, a positive definite matrix is a P -matrix, and every P -matrix is also a P0-matrix. For more properties about P -matrix and P0-matrix, please refer to[8]. The two concepts can be extended to nonlinear mappings.
Definition 2.2. Given a mapping F : Rn→ Rn, then F is a
(a) monotone function ifx − y, F (x) − F (y)0 for all x, y ∈ Rn; (b) P0-function if max1i n
xi=yi
(xi − yi)(Fi(x) − Fi(y))0 for all x, y ∈ Rnand x = y;
(c) P -function if max1i n(xi− yi)(Fi(x) − Fi(y)) > 0 for all x, y ∈ Rnand x = y;
(d) uniform P -function with modulus > 0 if max1i n(xi− yi)(Fi(x) − Fi(y)) x − y 2for all x, y ∈ Rn. From the above definitions, it is obvious that F is a P0-function if F is monotone, and the Jacobian matrix of every continuously differentiable P0-function is a P0-matrix. The following lemma states that the mapping Fεis a P -function if F is a P0-function.
Lemma 2.1 (Facchinei and Kanzow[9, Lemma 3.2]). For any ε > 0, let Fε : Rn → Rn be given by (10). If F is a P0-function, then the Jacobian matrices Fε (x) for all x ∈ Rn are P -matrices. In particular, the function Fε is a P -function.
Next, we review some favorable properties ofpwhere the proofs of Property 2.1 can be found in[2, Proposition 3.1]and[25, Lemmas 2.1 and 2.2]whereas the proof of Property 2.2 is given by[25, Lemma 3.1].
Property 2.1. Letp: R2→ R be defined as in (6). Then, the following results hold.
(a) pis an NCP-function.
(b) p is Lipschitz continuous with the Lipschitz constant L given by L =√
2+ 2(1/p−1/2) when 1 < p < 2 and L = 1 +√
2 when p 2.
(c) pis strongly semismooth.
(d) Given any point (a, b) ∈ R2, each element in the generalized gradient jp(a, b) has the representation ( − 1, − 1), where
=sgn(a) · |a|p−1
(a, b) p−1p and =sgn(b) · |b|p−1
(a, b) p−1p if (a, b) = (0, 0)
and otherwise (, ) ∈ R2denotes an arbitrary vector satisfying||p/(p−1)+ ||p/(p−1)1.
(e) If{(ak, bk)} ⊆ R2with ak → −∞, or bk → −∞, or ak→ +∞ and bk → +∞, then we have |p(ak, bk)| → +∞ for k → +∞.
Property 2.2. Letp : R2→ R be defined as in (6). Then, there exists two positive constants c1> 0 and c2> 0 such that c1| min{a, b}||p(a, b)|c2| min{a, b}|.
The semismooth property is very important from computational point of view. In particular, it plays a fundamental role in the superlinear convergence analysis of generalized Newton methods[26,28]. If the mapping G : Rn→ Rmis locally Lipschitz continuous, then G is almost everywhere differentiable by Rademacher’s Theorem (see[7]). In this case, the generalized JacobianjG(x) of G at x (in the Clarke sense) can be defined as the convex hull of the generalized JacobianjBG(x), where
jBG(x) := {V ∈ Rm×n|∃{xk} ⊆ DG: {xk} → x and G (xk) → V }
with DGdenoting the set of differentiable points of G. Assume that G : Rn→ Rmis locally Lipschitz continuous. G is called semismooth at x if G is directionally differentiable at x and for any V ∈ jG(x + h) and h → 0,
G(x + h) − G(x) − V h = o( h );
G is called strongly semismooth at x if G is semismooth at x and for any V ∈ jG(x + h) and h → 0,
G(x + h) − G(x) − V h = O( h 2); (16)
G is called a (strongly) semismooth function if it is (strongly) semismooth everywhere.
3. Properties ofHp(z) and Gp(z)
In this section, we will study the semismoothness of the mapping Hpand characterize its generalized Jacobian matrix at any point z. In particular, we also give a sufficient condition for the nonsingularity of all generalized Jacobians at a solution of (14). Then, we investigate some favorable properties of the merit function Gp(z) which are crucial to the convergence analysis of the regularized semismooth Newton algorithm described as in the next section.
Proposition 3.1. The mapping Hp : R+× Rn → Rn defined as in (13) is semismooth. Moreover, it is strongly semismooth if F is locally Lipschitz continuous.
Proof. Since a function is (strongly) semismooth if and only if its component functions are (strongly) semismooth, to prove that Hpis (strongly) semismooth we only need to prove that Hp,i, i = 1, 2, . . . , n + 1 are (strongly) semismooth.
Apparently, Hp,1is strongly semismooth by formula (16) since Hp,1(z) = ε. For Hp,i, i = 2, 3, . . . , n + 1, since p is strongly semismooth by Property 2.1 (c) and the composite of two (strongly) semismooth functions is (strongly) semismooth by[14, Theorem 19], we conclude that Hp,i, i = 2, 3, . . . , n + 1 are semismooth. If F is locally Lipschitz continuous, then Fεis strongly semismooth, and consequently, Hp,i, i = 2, 3, . . . , n + 1 are strongly semismooth.
We next give the estimation of the generalized Jacobian of Hpby Property 2.1 (d).
Proposition 3.2. For any z = (ε, x) ∈ R+× Rn, we have (jHp(z))T⊆
1 xTB(z)
0 (A(z) − I ) + (∇F (x) + εI )(B(z) − I )
, (17)
where A(z) and B(z) are possibly multi-valued n × n diagonal matrices with ith diagonal elements Aii(z) and Bii(z) given by
Aii(z) = sgn(xi) · |xi|p−1
(xi, Fε,i(x)) p−1p , Bii(z) =sgn(Fε,i(x)) · |Fε,i(x)|p−1 (xi, Fε,i(x)) p−1p
if (xi, Fε,i(x)) = (0, 0); and otherwise given by
Aii(z) = i, Bii(z) = i for any (i, i) such that |i|p/(p−1)+ |i|p/(p−1)1.
Proof. By the known rules on the evaluation of the generalized Jacobian (see[7, Proposition 2.6.2(e)]), we have jHp(z)T⊆ jHp,1(z) × jHp,2(z) × · · · × jHp,n+1(z),
where the right-hand side denotes a set of matrices whose ith column belongs to jHp,i(z), and Hp,iis the ith component function of Hp. Clearly,
jHp,1(z) =
1 0
∈ Rn+1.
For j = 2, 3, . . . , n + 1, letting i = j − 1 and applying Property 2.1 (d) yield jHp,j(z) =
sgn(xi) · |xi|p−1 (xi, Fε,i(x)) p−1p − 1
0 ei
+
xi
∇Fi(x) + εei
sgn(Fε,i(x)) · |Fε,i(x)|p−1 (xi, Fε,i(x)) p−1p − 1
if (xi, Fε,i(x)) = (0, 0); and otherwise jHp,j(z) = (i − 1)
0 ei
+
xi
∇Fi(x) + εei
(i− 1)
with|i|p/(p−1)+|i|p/(p−1)1, where eidenotes the vector whose ith element is zero and other elements are 1. From these equalities, the conclusion easily follows.
Now, exploiting the estimation ofjHp(z) given by (17), we may present a sufficient condition to guarantee the nonsingularity of all generalized Jacobians of Hpat a solution z∗of (14). This result is important for the superlinear (or quadratic) convergence of the semismooth Newton method (see[11]). Let z∗= (ε∗, x∗) ∈ R+× Rnbe a solution of (14). Clearly, ε∗= 0 and x∗is a solution of the NCP. For the sake of notation, let
I := {i ∈ {1, 2, . . . , n} | xi∗> 0, Fi(x∗) = 0}, J := {i ∈ {1, 2, . . . , n} | xi∗= 0, Fi(x∗) = 0}, K := {i ∈ {1, 2, . . . , n} | xi∗= 0, Fi(x∗) > 0}.
By rearrangement we assume that∇F (x∗) can be written as
∇F (x∗) =
∇FII(x∗) ∇FIJ(x∗) ∇FIK(x∗)
∇FJI(x∗) ∇FJJ(x∗) ∇FJK(x∗)
∇FKI(x∗) ∇FKJ(x∗) ∇FKK(x∗)
. (18)
The NCP is called R-regular at x∗ if ∇FII(x∗) is nonsingular and its Schur-complement in the matrix
∇FII(x∗)
∇FJI(x∗)
∇FIJ(x∗)
∇FJJ(x∗)
is a P -matrix.
Proposition 3.3. Suppose that z∗= (ε∗, x∗) ∈ R+× Rnbe a solution of (14) and the NCP is R-regular at x∗, then all V ∈ jHp(z∗) are nonsingular.
Proof. From Proposition 3.2, it is easy to see that for any V ∈ jHp(z∗)T, there exists a vector u(z∗) ∈ Rnand a matrix W (z∗) ∈ Rn×nsuch that
V =
1 u(z∗)T 0 W (z∗)
,
where
W (z∗) = (A(z∗) − I ) + (∇F (x∗) + ε∗I )(B(z∗) − I )
with A(z∗) and B(z∗) characterized as in Proposition 3.2. Therefore, proving that V is nonsingular is equivalent to arguing that W (z∗) is nonsingular. Using the expression of ∇F (x∗) in (18) and noting that ε∗= 0, we can rewrite W (z∗) in the partitioned form
W (z∗) =
−∇FII ∇FIJ(BJJ− IJJ) 0IK
−∇FJI ∇FIJ(BJJ− IJJ) + (AJJ− IJJ) 0JK
−∇FKI ∇FKJ(BJJ− IJJ) −IKK
,
where for convenience we dispense with the notations z∗ and x∗. The rest of the proof is identical to that of[11, Proposition 3.2].
In what follows, we concentrate on the properties of Gp. First, applying[6, Propositon 3.2 (c)]and Theorem 2.6.6 of[7], we immediately obtain the following conclusion.
Proposition 3.4. For any ε 0, the function p,ε defined by (12) is continuously differentiable everywhere, and consequently, Gpdefined as in (15) is continuously differentiable everywhere and∇Gp(z) = VTHp(z) for any V ∈ jHp(z).
Proposition 3.5. Suppose that F is a P0-function and ˆε, ˜ε are two given positive numbers such that ˆε < ˜ε. Then, the merit function Gpdefined as in (15) has the property:
k→+∞lim Gp(zk) = +∞
for any sequence{zk= (εk, xk)} such that εk∈ [ˆε, ˜ε] and xk → +∞.
Proof. We prove this by contradiction which is a standard and common technique. Suppose limk→+∞Gp(zk) = +∞.
Then from (15) and (12) it follows that there exists an unbounded sequence{xk} such that {p,εk(xk)} is bounded. Let J := {i ∈ {1, 2, . . . , n}|{xik} is unbounded}.
Since{xk} is unbounded, we have J = ∅. Without loss of generality, we assume that {|xjk|} → ∞ for any j ∈ J . Now, we define a bounded sequence by
yik :=
0 if i ∈ J, xik if i /∈ J.
From the definition of{yk} and F being a P0-function, we have 0 max
1i n xki=yki
(xik− yki)(Fi(xk) − Fi(yk))
= max
i∈J xik· (Fi(xk) − Fi(yk))
= xjk0· (Fj0(xk) − Fj0(yk)), (19)
where j0is one of the indices for which the max is attained. Since j0∈ J , we have that {|xjk0|} → +∞ as k → +∞. If xjk0 → −∞ as k → +∞, using Property 2.1(e) immediately yields that p(xjk0, Fεk,j0(xk)) → +∞. If xjk0 → +∞ as k → +∞, noting that Fj0(yk) is bounded by the continuity of Fj0, we have from (19) that Fj0(xk) does not tend to −∞, which in turn implies that{Fj0(xk)+εkxjk0} → +∞. From Property 2.1(e) where {xkj0} → +∞ and {Fj0(xk)+εkxjk0} → +∞,we also obtain that p(xjk0, Fεk,j0(xk)) → +∞. Thus, both cases yield p(xjk0, Fεk,j0(xk)) → +∞ which is a contradiction to the boundedness of{p,εk(xk)}. Consequently, we prove that limk→+∞Gp(zk) = +∞.
Remark 3.1. Proposition 3.5 implies thatp,εhas bounded level sets under the assumption of F being a P0-function.
However, from[6, Proposition 3.5], we know that a stronger condition (i.e., F being a uniform P -function) is needed to guarantee the level sets ofpto be bounded.
To close this section, we present two results which will be used to analyze the global convergence of the algorithm in the next section. The first result is extracted from Theorem 5.4 of[9], while the second result can be obtained by using Property 2.2 and following the same arguments as in[30, Proposition 2.2].
Proposition 3.6. Suppose that F is a P0-function and the solution set S∗of the NCP is nonempty and bounded. Suppose that {εk} and {xk} are two infinite sequences such that for each k 0, εk> 0, k0 satisfying limk→+∞εk = 0, limk→+∞k = 0. For each k 0, let xk ∈ Rn satisfy p(εk, xk) k. Then {xk} remains bounded and every accumulation point of{xk} is a solution of the NCP.
Proposition 3.7. Suppose that F is a monotone function and the solution set S∗of the NCP is nonempty. Suppose that {εk} and{xk} are two infinite sequences such that for each k 0, εk> 0, k0, kCεkand limk→+∞εk= 0, where C > 0 is a constant. For each k 0, let xk ∈ Rnsatisfy p(εk, xk) k. Suppose that x∗= arg minx∈S∗ x and F is Lipschitz continuous. Then{xk} remains bounded and every accumulation point of {xk} is a solution of the NCP.
4. Regularization semismooth Newton method
From the discussions of last section, we see that Hp(z) and Gp(z) for all p > 1 enjoy the same desirable properties.
Sun[30]used H2(z) and G2(z) to develop a regularization semismooth Newton method for the NCP. In this section, we will develop a regularization semismooth Newton algorithm by any Hp(z) and Gp(z) with p > 1. This algorithm is guaranteed to solve P0-complementarity problems due to Proposition 3.5.
Now we are ready to describe this specific algorithm. We adopt almost the same notations used in[30]. Choose
¯ε ∈ (0, +∞) and ∈ (0, 1) such that ¯ε < 1. Let t ∈ [1/2, 1] and ¯z := (¯ε, 0) ∈ R++×Rn. Define : R+×Rn→ R+ by
(z) := min{1, Gp(z)t}. (20)
We also denote
:= {z = (ε, x) ∈ R+× Rn|ε(z)¯ε}. (21)
Note that(z) for any z ∈ R+× Rnby (20). Hence, (¯ε, x) ∈ for any x ∈ Rn. In addition, by the definition of (z), it is easily shown the following relation holds.
Proposition 4.1. Let Hpand be defined as in (13) and (20), respectively. Then, Hp(z) = 0 ⇐⇒ (z) = 0 ⇐⇒ Hp(z) = (z)¯z.
Algorithm 4.1 (The Regularization Newton Algorithm).
(Step 0) Given any p > 1 and choose constants ∈ (0, 1), t ∈ [1/2, 1] and ∈ (1, 1/2). Let ε0:= ¯ε and x0∈ Rnbe an arbitrary point. Set k := 0.
(Step 1) If Hp(zk) = 0, then stop. Otherwise, let k:= (zk) = min{1, Gp(zk)t}.
(Step 2) Choose Vk∈ jHp(zk) and compute zk= (εk, xk) ∈ R × Rnby
Hp(zk) + Vkzk= k¯z. (22)
(Step 2) Let lkbe the smallest nonnegative integer l such that
Gp(zk+ lzk)[1 − 2 (1 − ¯ε) l]Gp(zk). (23)
(Step 2) Set zk+1 := zk+ lkzk. (Step 4) Set k := k + 1 and go to Step 1.
From Proposition 3.2, we know that for any V ∈ jHp(z) with z=(ε, x) ∈ R++×Rn, there exists a W =(u(z) W (z)) ∈ jp(z) with u(z) ∈ Rnand W (z) ∈ Rn×nsuch that
V =
1 0
u(z) W (z)
. (24)
Suppose that F is a P0-function. Then by Lemma 2.1 Fε (x) is a P -matrix. Hence, for any x ∈ Rnand ε > 0, W (z) is nonsingular by the proof of Proposition 2 of[20]. It thus follows that all V ∈ jHp(z) with z = (ε, x) ∈ R++× Rn are nonsingular. Therefore, the Newton step in (22) is well-defined, and moreover, from (22), for any k 0 and εk> 0, there exists a Wk∈ jp(zk) such that
(∇p(zk))Tzk= p(zk)TWkzk= −p(zk)Tp(zk) = −2p(zk). (25) Using the equality and Proposition 4.1, we next show that Algorithm 4.1 is well-defined.
Proposition 4.2. Suppose that F is a P0-function and zk= (εk, xk) ∈ R++× Rnfor k 0. Then zk+1 ∈ R++× Rn and Algorithm 4.1 is well-defined.
Proof. Since εk k= (zk) > 0. From the first component in the relation (22) in Algorithm 4.1, we have
εk+ εk= k¯ε ⇒ εk= −εk+ k¯ε. (26)
Then, for any ∈ [0, 1], there has
εk+ εk= (1 − )εk+ k¯ε > 0. (27)
Thus, combining the fact that(z) Gp(z)1/2with (22) and (27) yields that (εk+ εk)2= [(1 − )εk+ k¯ε]2
= (1 − )2(εk)2+ 2(1 − )kεk¯ε + 22k¯ε2
(1 − )2(εk)2+ 2kεk¯ε + O(2)
(1 − )2(εk)2+ 2 Gp(zk)1/2 Hp(zk) ¯ε + O(2)
= (1 − 2)(εk)2+ 2√
2 ¯εGp(zk) + O(2). (28)
Now, we define
() := p(zk+ zk) − p(zk) − (∇p(zk))Tzk.
Sincepis continuously differentiable at any zk ∈ R++× Rnby Proposition 3.4, we obtain() = o(). On the other hand, from (22) and (25) it follows that
1
2 p(zk+ zk) 2= p(zk+ zk)
= p(zk) + (∇p(zk))Tzk+ ()
= p(zk) − 2p(zk) + o()
= (1 − 2)p(zk) + o() (29)
for any ∈ [0, 1]. Therefore, using Eqs. (28) and (29), we obtain Gp(zk+ zk) =12 Hp(zk+ zk) 2
=12(εk+ εk)2+12 p(zk+ zk) 2
12(1 − 2)(εk)2+√
2 ¯εGp(zk) + (1 − 2)p(zk) + o()
(1 − 2)Gp(zk) + 2 ¯εGp(zk) + o()
= [1 − 2(1 − ¯ε)]Gp(zk) + o() (30)
for any ∈ [0, 1]. The inequality (30) implies that there exists ¯ ∈ (0, 1] such that Gp(zk+ zk)[1 − 2 (1 − ¯ε)]Gp(zk) ∀ ∈ [0, ¯],
which indicates that Algorithm 4.1 is well-defined.
Proposition 4.3. Suppose that F is a P0-function. For each k 0, if εk> 0 and zk ∈ , then for any ∈ [0, 1] such that
Gp(zk+ zk)[1 − 2 (1 − ¯ε)]Gp(zk), (31)
there holds that zk+ zk∈ .
Proof. We prove this proposition by considering the following two cases:
Case (i): Gp(zk) > 1. Then k= . From zk ∈ and (z) = min{1, Gp(z)t} for any z ∈ R+× Rn, it follows that for any ∈ [0, 1],
(εk+ εk) − (zk+ zk)¯ε (1 − )εk+ k¯ε − ¯ε
(1 − )k¯ε + k¯ε − ¯ε
= 0. (32)
Case (ii): Gp(zk)1. Then, for any ∈ [0, 1] satisfying (31), we have
Gp(zk+ zk)[1 − 2 (1 − ¯ε)]Gp(zk)1. (33)
Therefore, for any ∈ [0, 1] satisfying (31), (zk+ zk) = Gp(zk+ zk)t.
Using the fact that zk ∈ and the first inequality in (33), we then obtain that for any ∈ [0, 1] satisfying (31), (εk+ εk) − (zk+ zk)¯ε (1 − )εk+ k¯ε − Gp(zk+ zk)t¯ε
(1 − )k¯ε + k¯ε − [1 − 2 (1 − ¯ε)]tGp(zk)t¯ε
= k¯ε − [1 − 2 (1 − ¯ε)]tGp(zk)t¯ε
= Gp(zk)t¯ε − [1 − 2 (1 − ¯ε)]tGp(zk)t¯ε
= {1 − [1 − 2 (1 − ¯ε)]t}Gp(zk)t¯ε
0. (34)
Combining (32) and (34) immediately yields the desired result.
Proposition 4.4. Suppose that F is a P0-function. Then Algorithm 4.1 generates an infinite sequence{zk} with zk ∈ for all k and
0 < εk+1εk ¯ε for all k. (35)
Proof. Since z0= (¯ε, x0) ∈ , the first part of the conclusions follows by repeatedly resorting to Propositions 4.2 and 4.3. We next concentrate on the proof of (35). First, ε0= ¯ε > 0. From the design of Algorithm 4.1 and the fact that (z) = min{1, Gp(z)t} for any z ∈ R+× Rn, it then follows that
ε1= (1 − l0)ε0+ l0(z0)¯ε (1 − l0)¯ε + l0 ¯ε ¯ε.
Hence (35) holds for k = 0. Suppose that (35) holds for k = i − 1. We next prove that (35) holds for k = i. From the design of Algorithm 4.1, we have
εi+1= (1 − li)εi+ li(zi)¯ε.
Noting that εi(zi)¯ε since zi ∈ , we then obtain εi+1(1 − li)εi + liεi= εi
and
εi+1(1 − li)(zi)¯ε + li(zi)¯ε = (zi)¯ε > 0.
Therefore, (35) holds for k = i. We complete the proof.
Now, using Propositions 3.5–3.7 and Proposition 4.4 and following the same arguments as in[30], we obtain the following global convergence results of Algorithm 4.1.
Theorem 4.1. Suppose that F is a P0-function and the solution set S∗of the NCP is nonempty and bounded. Then the infinite sequence{zk} generated by Algorithm 4.1 is bounded and any accumulation point of {zk} is a solution of H (z) = 0.
Theorem 4.2. Suppose that F is a monotone function and in Algorithm 4.1 the parameter t =12. Then if the iteration sequence{zk} is bounded, then the solution set S∗of the NCP is nonempty. Conversely, if the solution set S∗of the NCP is nonempty and F is Lipschitz continuous, then the infinite sequence {zk} generated by Algorithm 4.1 is bounded and any accumulation point of{zk} is a solution of H(z) = 0.
In addition, using Proposition 3.1 and similar proof as for[30, Theorem 5.1], we obtain the following local superlinear (quadratic) convergence results of Algorithm 4.1.
Theorem 4.3. Suppose that F is a P0-function and the solution set S∗of the NCP is nonempty and bounded. Suppose that z∗ := (ε∗, x∗) ∈ R × Rnis an accumulation point of the infinite sequence{zk} generated by Algorithm 4.1 and all V ∈ jHp(z∗) are nonsingular. Then the whole sequence {zk} converges to z∗with
zk+1− z∗ = o( zk− z∗ ), εk+1= o(εk).
Furthermore, if F is locally Lipschitz continuous around x∗, then zk+1− z∗ = O( zk− z∗ 2), εk+1= O(εk)2.
Moreover, from Proposition 3.3, all the conclusions of Theorem 4.3 hold if the assumption that all V ∈ jHp(z∗) are nonsingular is replaced by that the NCP is R-regular at x∗.
5. Numerical experiments
We implemented Algorithm 4.1 by our codes in MATLAB 6.5 for almost all test problems except the unavailable
“pvg105” and “scarfbnum” with the starting points in MCPLIB[1]. All numerical experiments were done at a PC with CPU of 2.8 GHz and RAM of 512 MB. Throughout the experiments, unless otherwise stated, we adopted the following parameters for Algorithm 4.1:
= 0.5, t = 1/2, = 10−4, = 0.5, ¯ε = 0.1.
0 20 40 60 80 100 102.02
102.03 102.04 102.05 102.06 102.07
value of p
the number of iteration
Number of iteration v.s. different p
Fig. 1. The number of iterations vs. value of p for Example “bertsekas(3)”.
We terminated the iteration if one of the following conditions was satisfied:
(1) Hp(zk) 1and min{xk, F (xk)} 2; (2) the step lengthk= lk is less thanmin. (3) the number of iteration exceeds kmax.
Among others, in our implementation the termination parameters were chosen as follows:
1= 10−10, 2= 10−6, min= 10−25 and kmax= 1000.
During the experiments, we incorporated some strategies to improve the numerical behavior of Algorithm 4.1 to some extent. These strategies are well-accepted and used in basically all suitable implementations of complementarity solvers.
The first modification is in the line search step. We replaced the standard (monotone) Armijo-rule by nonmonotone line search described in[35]to seek a suitable steplength, i.e., we computed the smallest nonnegative integer l such that
Gp(zk+ ldk)Wk− 2 (1 − ¯ε) lGp(zk) for all k 0, where Wkis given by
Wk= (k−1Qk−1Wk−1+ Gp(zk))/Qk
with
Qk= k−1Qk−1+ 1.
In our implementation, we usedW−1= Gp(z0), Q−1= 1, −1= 0.85 and k ≡ 0.85.
The second modification is necessary since the mapping F is often not defined outside the positive orthant whereas our algorithm assumes that F can be evaluated on the whole space Rn. Hence, in order to avoid possible domain violations, we employed a simple backtracking strategy: Given an iterate zk = (εk, xk) ∈ R++× Rn and a search direction dk ∈ Rn+1, we first compute the exponent jk := min{0, 1, 2, . . . , } such that
F (xk+ jkdk(2 : n + 1))
