Numerical results for a globalized active-set Newton method for mixed complementarity problems

Numerical results for a globalized active-set Newton method for mixed complementarity problems

A.N. Daryina^I; A.F. Izmailov^II; M.V. Solodov^III

^IRussian Peoples Friendship University, Miklukho-Maklaya Str. 6, 117198 Moscow, Russia, E-mail: [e-mail]

^IIMoscow State University, Faculty of Computational Mathematics and Cybernetics, Department of Operations Research, Leninskiye Gori, GSP-2, 119992 Moscow, Russia, E-mail: [e-mail]

^IIIInstituto de Matemática Pura e Aplicada, Estrada Dona Castorina 110, Jardim Botânico 22460-320 Rio de Janeiro, RJ, Brazil E-mail: [e-mail]

ABSTRACT

We discuss a globalization scheme for a class of active-set Newton methods for solving the mixed complementarity problem (MCP), which was proposed by the authors in [3]. The attractive features of the local phase of the method are that it requires solving only one system of linear equations per iteration, yet the local superlinear convergence is guaranteed under extremely mild assumptions, in particular weaker than the property of semistability of an MCP solution. Thus the local superlinear convergence conditions of the method are weaker than conditions required for the semismooth (generalized) Newton methods and also weaker than convergence conditions of the linearization (Josephy-Newton) method. Numerical experiments on some test problems are presented, including results on the MCPLIB collection for the globalized version.

Mathematical subject classification: 90C30, 65K05.

Keywords: mixed complementarity problem, semistability, 2-regularity, weak regularity, error bound, Newton method.

1 Introduction

The mixed complementarity problem (MCP) [9] is the variational inequality on a generalized box, that is

where F: ⁿ® ⁿ and

B = {x Î ⁿ |l_i < x_i < u_i, i = 1, ¼, n},

l_i Î È{-¥}, u_i Î È{+¥}, l_i < u_i for all i = 1, ¼, n. Equivalently, it can be stated as

find x Î B such that F_i(x)

As is well known, many important problems can be cast in the format of MCP [11, 9]. As a special case of MCP, we mention the nonlinear complementarity problem (NCP), which corresponds to setting l_i = 0, u_i = +¥, i = 1, ¼, n. The systems of nonlinear equations are obtained by choosing l_i = -¥, u_i = +¥, i = 1, ¼, n. Another important example is the primal-dual Karush-Kuhn-Tucker (KKT) optimality system: find z Î ^p and µ Î ^m such that

where g: ^p® ^p and G: ^p® ^m. The KKT system (1.2) can be written as an MCP if we set n = p + m and

F(x) = , x = (z, µ) Î ^p × ^m,

l_i = -¥, i = 1, ¼, p, l_i = 0, i = p + 1, ¼, n, u_i = +¥, i = 1, ¼, n.

Under well-known assumptions, (1.2) represents the first-order primal-dual necessary conditions characterizing solutions in variational inequality or inequality-constrained optimization problems.

This paper describes a globalization scheme for a local Newton-type method for solving MCP, which was proposed by the authors in [3]. Also numerical experiments will be reported on problems from the MCPLIB collection [4], together with results on some additional test examples.

The algorithm of [3] belongs to the class of active-set methods, which isrelated to the idea of identifying active constraints in (inequality-)constrained optimization. The importance of identifying active constraints in optimization is well recognized, since it allows one to locally reduce the problem to an equality-constrained one, the latter being structurally significantly easier. For some identification techniques we refer the reader to [8, 17, 6] for the case of nonlinear constraints, and to [5, 7] for linear (complementarity) problems. The linear case, of course, is very special: correct identification is easier and can often be obtained (relatively) far from a solution. In the nonlinear case, correct identification is in general quite local. That said, optimization algorithms usually try to identify the correct set ''early'' (even if the identification cannot be provably correct). The intended use of identification in our method is somewhat different. We view it as a strictly local phase and see no reason of even trying it when far from a solution (i.e., outside of the region where forcing quadratic convergence can be expected). For this purpose, our method has some checks that prevent switching to the local active-set phase, unless there are reasons to believe that we are really close to a solution.

A few words about our notation. Given a finite set I, |I| stands for its cardinality. By (m, n) we denote the space of m × n matrices with real entries. By E we shall denote the identity matrix whose dimension would be always clear from the context. For x Î ⁿ and an index set I Ì {1, ¼, n} , x_I stands for the vector with components x_i, i Î I. For a linear operator L, im L is its range (image space), and ker L is its kernel (null space). For a directionally differentiable mapping f: ⁿ® ^m, by f¢(x; d) we denote the usual directional derivative of f at x Î ⁿ in the direction d Î ⁿ. If {z^k} is a sequence in ^p and {t_k} is a sequence in such that t_k ® 0+ as k ® ¥, by z^k = o(t_k) we mean that lim_{k ® ¥}||z^k||/t_k = 0.

2 The local active-set method

In the context of MCP, active-set strategy corresponds to identifying the sets of indices

where is some solution of MCP. If the specified sets can be correctly identified using information available at a point x close enough to the solution , then locally MCP can be reduced to a system of nonlinear equations (which is structurally a much simpler problem to solve). In the sequel, we shall also use the following partitioning of the set of active indices:

The analog of the strict complementarity condition in NCP (or KKT) corresponds, in the setting of MCP, to saying that A₀ = Æ. Under this assumption, locally MCP trivially reduces to a system of nonlinear equations, which simplifies the local structure of MCP significantly. The condition of strict complementarity, however, is restrictive and will not be assumed.

Let y: × ® be a complementarity function, i.e., a function such that

Assume that y satisfies the following additional assumptions:

Then MCP can be equivalently reformulated as a system of nonlinear equations

where

For some special choices of complementarity functions, this reformulation is well-known (see, e.g., [1, 10]).

Complementarity functions to be used in the sequel are the natural residual

y_NR(a, b) = min{a, b},

the Fischer-Burmeister function

and

y_S(a, b) = 2ab - (min{0, a + b})² ,

where S stands for ''smooth''. Note that each of the three functions satisfies the conditions stated above.

The corresponding reformulations of MCP would be denoted by Y_NR, Y_FB and Y_S, respectively.

The identification of relevant active index sets would be based on error bound analysis, which we describe next.

Definition 2.1.

Let Y:

ⁿ ® ⁿ be differentiable in a neighbourhood of Î ⁿ and Y¢ : ⁿ® (n, n) be directionally differentiable at . Then Y is 2-regular at if

T = {0} ,

where

with P being the orthogonal projector onto (im Y¢())^{^}.

The above is a special case of 2-regularity of a nonlinear mapping [13, 12] corresponding to the case when the mapping acts from some space into itself.

Theorem 2.1 [3, Theorem 2.2]. Let F :

ⁿ® ⁿ be sufficiently smooth near a point Î ⁿ , which is a solution of MCP.

The mapping Y_S is 2-regular at

if, and only if, there exist a neighborhood U of

and a constant M > 0 such that

Adjusting M and U, if necessary, the error bound (2.3) can be simplified into the following relation (less accurate, but usually easier to use):

The assumption of 2-regularity above is extremely mild. In particular, it is strictly weaker than various alternatives, such as semistability [2] of the MCP solution (which is already a very mild condition, weaker than strong regularity or quasi-regularity). Specifically, we have the following.

Proposition 2.1 [3, Proposition 2.3]. Semistability of a solution

of MCP implies 2-regularity of Y_S at

(equivalently, error bound (2.3)), but not vice versa.

As a consequence, the error bound given above can hold even when MCP reformulations based on Y_NR and Y_FB fail to provide a bound (the latter bounds can be shown to be equivalent to semistability).

The following technique for identifying the relevant index sets is based on the ideas of [6], see also [9, Ch. 6.7]. Define the identification function

where Î (0, 1) and > 0 are fixed numbers (the choice of and does not affect theoretical analysis; in our numerical experiments reported in Section 3, we use = 0.9 and = -1/log, as suggested in [6]). For any x Î ⁿ, define further the index sets

Proposition 2.2 [3, Proposition 3.1]. If Y_S is 2-regular at a solution

of MCP (equivalently, the error bound (2.3) holds), then for any x Î ⁿ sufficiently close to

, it holds that

Observe that in the implementation of the identification procedure, the following obvious relations can also be taken into account:

Once the index sets are correctly identified, we have the following relations which are guaranteed to be satisfied at a solution of MCP:

For simplicity of notation, suppose that the components of x Î ⁿ are ordered in such a way that x = . Then MCP locally reduces to the following system of nonlinear equations:

Observe that in the absence of strict complementarity (when A₀¹ Æ, i.e.,|A| > |A₊|), the system is over-determined (the number of equations is larger than the number of unknowns). This opens up a number of options. Of course, one can just solve the system by the Gauss-Newton method (GNM). This possibility will be considered. However, we prefer not to limit ourselves to GNM for the following reason: the Gauss-Newton approach can destroy structure present in F_A (for example, sparsity or the primal-dual structure in the case of KKT).

Our proposal is to consider the following system of nonlinear equations:

where

with C : ^|A+|® (|A₊|, A) being a smooth mapping (possibly constant). Clearly, is a solution of (2.13) for any choice of C. The Jacobian of (2.13) at this solution is given by

where we have taken into account that F_A() = 0. Thus can be found by applying Newton-type methods to (2.13) whenever the matrix in (2.13) is nonsingular.

Note that GNM for (2.12) would essentially correspond to choosing in (2.13)

and applying to the resulting system an approximate version of the pure Newton method. Indeed, with the notation of (2.15), the Gauss-Newton iteration for (2.12) has the form

Observe that the above formula is just an approximation of the standard Newton iteration for (2.13), where the Jacobian is replaced by

Due to (2.14), this change preserves the superlinear convergence of the pure Newton iteration for (2.13). Note finally that with the choice of (2.15), we have

This immediately motivates the following definition.

Definition 2.2. A solution of MCP is referred to as weakly regular if

Weak regularity is implied by semistability, but not vice versa. Moreover, 2-regularity of Y_S at and weak regularity, when combined, are still a strictly weaker condition than semistability.

Proposition 2.3 [3, Proposition 3.3]. Let

be a solution of MCP. Then semistability of

implies weak regularity of

, but not vice versa.

Proposition 2.4 [3, Proposition 3.4]. Let

be a solution of MCP. Then semistability of

implies the combination of 2-regularity of Y_S at

and weak regularity of

, but not vice versa.

We have thus a local algorithm with superlinear convergence under assumptions weaker than semistability of the MCP solution. Specifically, we have the following.

Theorem 2.2 [3, Theorem 3.5]. Let F :ⁿ ®

ⁿ be sufficiently smooth near a point Î ⁿ , which is a solution of MCP. Suppose that this solution is weakly regular and Y_S is 2-regular at

.

For any x⁰Î ⁿ sufficiently close to

, if the index sets A = A(x⁰), A₊ = A₊(x⁰), A_0l= A_0l(x⁰), A_0u= A_0u(x⁰), N_l= N_l(x⁰) and N_u= N_u(x⁰) are defined according to (2.5)-(2.9), then GNM applied to the system (2.12)(with

as a starting point) is well-defined and superlinearly convergent to .

As already mentioned above, it sometimes can be useful to choose the mapping C differently from the Gauss-Newton option of (2.15). We might want to take C(·) = C Î (|A₊|, |A|), a fixed matrix, in order to preserve in the matrix

the structure (primal-dual, sparsity, etc.) of the matrix

This motivates the following considerations.

Proposition 2.5 [3, Proposition 3.6]. Suppose that a solution

of MCP is weakly regular.

Then the set of matrices C Î (|A₊|,|A|) such that () is nonsingular is open and dense in (|A₊|,|A|).

In the situation of weak regularity of the solution, Proposition 2.5 justifies choosing C in any desirable way, as the chance that the resulting system would be degenerate is negligible (the set of matrices for which this would happen is of the Lebesgue measure zero). Of course, one should make reasonable choices. For example, it should hold that rankC = |A₊|.

3 Globalization issues and numerical experiments

In this section we report on our numerical experience based on the MCPLIB test problems collection (the newer version of [4]), and on some additional small examples, designed to highlight the case where various standard regularity conditions do not hold, and thus the semismooth (generalized) Newton based methods (SNM) may have trouble or converge slowly. This is precisely the case where the switch to our local algorithm can be particularly useful.

To perform numerical experiments, we had to implement our method as a final stage of some globally convergent scheme. It seems difficult to suggest a globally convergent scheme directly related to the structure of our local algorithm. In some sense, this is a disadvantage. But on the other hand, our local approach can be combined with any globally convergent algorithm satisfying some requirements (see below). In fact, the way we suggest to use our active-set method is precisely for improving the local convergence properties of standard algorithms (say, when they run into difficulties because of the lack of regularity of a solution). Our numerical experiments, reported below, indicate that the resulting strategy fulfills the stated objective.

We now describe the basic hybrid globalization scheme and its convergence properties. The scheme employs the linear decrease test for some (computable) merit function j : ⁿ ® ₊. Possible choices of a global algorithm and j are restricted by the following assumptions:

(A1) j(x) = 0 if, and only if, x Î ⁿ is a solution of MCP.

(A2) j decreases monotonically along the trajectories of the global algorithm.

(A3) j decreases superlinearly along the trajectories of our local method near "qualified" solutions, i.e., solutions satisfying the assumptions of the local convergence theorem for our method (Theorem 2.2).

Our choice of a global algorithm and a merit function, as well as the full method, would be stated later. We first give a simplified general discussion.

At each iteration k of the global algorithm, we first identify the relevant index sets. If there is a reason to believe that the identification is correct (e.g., if it did not change when compared to the previous iteration), then we compute the trial point by the step of GNM applied to (2.12) at , and set the remaining components of Î ⁿ equal to the corresponding components of and . If is well-defined and

then we set x^k+1 = , and proceed with the next iteration. Otherwise, we compute x^k+1 by the step of the global algorithm, and proceed with the next iteration.

If the linear decrease test (3.1) is satisfied for a finite number of iterations only, the hybrid globalization scheme behaves essentially as the global algorithm, and the global convergence properties of the latter remain valid. Otherwise, taking into account assumption (A2), we conclude that

j(x^k+1) ® 0 as k ® ¥,

and in particular, according to assumption (A1), each accumulation point of the sequence {x^k} is a solution of MCP. This characterizes the global convergence properties of the hybrid scheme.

Suppose now that a ''qualified'' solution of MCP is an accumulation point of {x^k} . Then, by assumption (A3), the linear decrease test (3.1) is satisfied for all iteration indices k sufficiently large. Furthermore, it is easy to see that the hybrid scheme eventually switches to our local algorithm (with the ''starting point'' sufficiently close to ). Hence, by Theorem 2.2, the entire sequence {x^k} converges to superlinearly.

In our numerical experiments we adopt the following choices. As a global algorithm, we use the linesearch-SNM for (2.1) with Y = Y_FB. More precisely, we essentially follow the implementation suggested in [15] (''General LineSearch Algorithm''), with all the parameter values adopted there. Furthermore, we take j = j_FB, where

We emphasize that this is just one example of the various appropriate choices. But this particular algorithm seems reasonable for our purposes. Linesearch SNM is known to be quite efficient, and at the same time, it is easy to implement in its basic form. We note that we do not use any enhancements, such as crashing and nonmonotone linesearch (see, e.g., [16]). The reason is that these are intended to improve global behavior of the algorithm, while we are concerned with local behavior. Thus a simple implementation of the global scheme is sufficient for our purposes, as our principal conclusions refer to the local convergence properties.

Assumptions (A1) and (A2) are evidently satisfied for the adopted choices. In order to guarantee (A3), we need to assume the error bound

Recall that this bound is equivalent to semistability. Under this assumption, for x⁰Î ⁿ close enough to a ''qualified'' solution , and the corresponding trajectory {x^k} of the local method, we have

where the second equality is by the Lipschitz-continuity of Y_FB near , the third is by Theorem 2.2, and the fourth is by (3.2). Thus, assumption (A3) is satisfied.

We proceed with the formal statement of the algorithm. Variable ''Alg'' below is used to select between the two variants of the algorithm that would be compared to each other. Alg = 1 means that the possibility of switching to our step is blocked, and so the algorithm works as the usual linesearch SNM/FB from [15]. Alg = 2 corresponds to the proposed active-set strategy. Note however that switching to an active-set step is forbidden on the first iteration, and also in the case when the identified index sets differ from the corresponding index sets at the previous iteration. This is done in order to prevent the algorithm from switching to the active-set strategy too early, when the sets are not yet stabilized and are likely to give incorrect identification.

Algorithm 3.1

Preliminary step. Set Alg = 1 or 2. Fix q, e, t Î (0, 1), d, g > 0. Set k = 0 and choose x⁰Î ⁿ.

Initialization step. If Alg = 1 or k = 0, go to SNM/FB step. Otherwise define the index sets A = A(x^k), N_l = N_l(x^k), N_u = N_u(x^k), A₊ = A₊(x^k), A_0l = A_0l(x^k), A_0u = A_0u(x^k) according to (2.5)-(2.9). If at least one of these sets does not coincide with its counterpart computed at the previous iteration, go to SNM/FB step.

GNM/active-set step. Compute Î ⁿ as follows:

If this point is well-defined and

set x^k+1 = , adjust k by 1, and go to Initialization step.

SNM/FB step. Compute L_k Î ¶_BY_FB(x^k) and

If this point is well-defined and (3.4) holds, set x^k+1 = , adjust k by 1, and go to Initialization step.

If is well-defined but (3.4) does not hold, set = - x^k. If

set d^k = and go to Linesearch step.

Gradient step. Set d^k = -(x^k).

Linesearch step. Compute the stepsize parameter a_k according to the Armijo rule: a_k = t^s, where s is the smallest nonnegative integer satisfying

j_FB(x^k + t^sd^k) < j_FB(x^k) + et^s áj¢_FB(x^k), d^kñ.

Set x^k+1 = x^k + a_kd^k, adjust k by 1, and go to Initialization step.

Note that (3.3) differs slightly from the pure Gauss-Newton iteration given by (2.16) with C(x_A+) defined in (2.15). The modification is made in order to reduce the number of evaluations of the Jacobian of F: the described algorithm requires exactly one such evaluation per iteration, whether Alg = 1 or 2. This modification does not affect the rate of convergence. To see this, note that (2.16) can be written in the form

Comparing this iteration with (3.3), observe that

provided the index sets are correctly identified for a given solution .

To compute L_kÎ ¶_BY_FB(x^k) at SNM/FB step, we used the procedure suggested in [1] (with z_i = 1 " i = 1, ¼, n, in the notation of [1], and the ''computer zero'' parameter set to 10^-10). The well-known formula

is used to compute the gradient of j_FB.

In the numerical experiments reported below, we used the following set of parameters: q = 0.9, e = 10^-4, t = 0.5, d = 2.1, g = 10^-9. The stopping criterion is

The cases when an algorithm did not terminate according to this criterion after 500 iterations are referred to as failures.

The algorithm was implemented in Matlab, making use of the standard option for treating sparse matrices. The Table 1 below reports first the total number of iterations for the two algorithms, and then the number of times the GNM/active-set step has been used. Next to the latter number, in the brackets, stands the number of active-set iterations at the ''tail'' of the process (right before convergence had been declared). Next, the total number of SNM/FB steps of the globalization scheme is reported for the two algorithms. In the brackets is the number of complete SNM/FB steps, i.e., those accepted with the unite stepsize in the linesearch step or by the linear decrease test without linesearch. Finally, the number of gradient steps and the number of evaluations of F are reported (recall that for both algorithms, the Jacobian of F is evaluated once per iteration). Failures are marked by ''-''.

Both variants of the algorithm failed for the following 10 problems: billups, bishop, colvdual, ehl_k40, ehl_k60, ehl_k80, forcebsm, forcedsa, lincont, simple-ex. Recall that by failure we mean that (3.5) was not satisfied after 500 iterations. In particular, for billups both variants of the algorithm converged to a local minimizer of the merit function, which is not an MCP solution. But for this problem, this is typical for (most) MCP algorithms. We note that the problems mentioned above are considered among the difficult ones in MCP literature. Since we were able to solve essentially all the other problems, this means that our implementation of the global scheme, though simple, is sufficiently robust. The rest of this discussion focuses on the cases for which at least one algorithm did not fail.

For 3 problems, namely duopoly, shubik, and ne-hard, switching to our step at some early stage prevented failure (although for ne-hard, the accuracy obtained by the algorithm without switching was of order 10^-9, i.e., almost satisfying the stopping test (3.5)). The opposite situation was observed for games only, though the obtained accuracy in the latter case was also of order 10^-9. We do not have an explanation for the behavior on those 4 problems and regard it as (probably) ''accidental''. Moreover, the possible negative global effect of switching to the active-set step too early can be avoided by the following simple trick. When an active-set step is accepted, we can store the previous iterate as a back-up, and restart the algorithm from that point if an active-set step is rejected on some subsequent iteration (which indicates that the switch occurred too early).

For 20 test problems the active-set step was never accepted, but trying it never harmed drastically. The number of evaluations of F for the algorithm without the switching option is typically not much less than for the complete variant, especially when the number of iterations is relatively large. We consider this extra work as a price to pay for safeguarding better local convergence properties in situations where SNM/FB runs into difficulties (see below).

The GNM/active-set step was ever accepted ''improperly'', i.e., far from a solution and with (probably) incorrect identification, for 5 test problems only. As mentioned above, for duopoly and shubik this actually prevented failure, but for games this caused it. For freebert this resulted in some extra linesearch steps. This points to the (obvious) fact that switching to the active-set strategy too early should be avoided.

For 19 test problems our step was used ''properly'' (on the final stage of the process). For badfree and handskoop, this was evidently rather advantageous. This has to do, of course, with the lack of (BD-)regularity, which affects the SNM/FB algorithm. In the other cases (apparently either BD-regular or not regular even in our sense), the conclusions are overall similar to the situation when our step was never used. In those cases, we do not win or lose much. On the one hand, we typically have to pay the price of a few more evaluations of F (but not always). On the other hand, note that the dimension of the linear system of equations which gives the GNM/active-set step is smaller than the dimension of the system to compute the SNM/FB step. This can be significant for large-scale problems (if the GNM/active-set step is accepted then no further computation is needed at this iteration). Finally, the ''proper'' use of our step never increased (and sometimes decreased) the overall number of iterations.

General conclusions are as follows. The option of switching to our step never harms too much, though we certainly have to pay some extra price for computing it at some iterations (at those where the index sets do not change), even if this step is eventually rejected. But this is consistent with the main goal of the presented approach. We emphasize that the goal is not in improving the linesearch SNM/FB (or any other algorithm) when it works efficiently. We try not to harm/interfere too much in those cases, while extending the algorithm to the (irregular) cases when SNM/FB does not work well.

We also point out that the efficiency of the identification procedure is of crucial importance for the methods presented in this paper. According to our numerical experience, the switching to our local method usually occurs exactly one iteration after the correct identification is obtained. Additional tuning of the identification procedure (i.e., using scaling or different identification functions) might improve the performance significantly. Actually, identification techniques other than the one described above could also be tried. Also, note that we have not been using any heuristic considerations to decide whether or not the active-set step should be computed. Developing such heuristics can certainly save some computational work as well. For example, even if the active sets have not changed from one iteration to the next, we may decide not to compute the active-set step if the residual is relatively large (i.e., we are still far from a solution). Other important issues are feasible versions of the method, and different globalization schemes which would better fit the structure of the method. This will be the subjects of future research.

To conclude, we illustrate some possible scenarios of the purely local behavior of SNM/FB and GNM/AS by applying them to the following small test problems with violated BD-regularity for Y_FB. We are talking here about the basic SNM/FB and GNM/AS iterations, without any modifications and tricks concerned with globalization. In particular, GNM/AS algorithm is specified precisely by (2.15), (2.16), with the index sets identified at the starting point (thus, this is the algorithm whose local properties are characterized by Theorem 2.2).

The first problem is a slight modification of [3, Example 3.1].

Example 3.1. Let n = 2, l_i = 0, u_i = +¥, i = 1, 2, and let F(x) = ((x₁ - 1)², x₁ + x₂ + -1). The point = (1, 0) is the solution of this NCP. Semistability (and hence, BD-regularity for Y_NR and Y_FB) is violated here, while 2-regularity and weak regularity hold. The starting point is x⁰ = (1.5, -0.5), with ||x⁰ - || » 7.1e-01, ||Y_FB(x⁰)|| » 8.2e-01, det L₀» 1.4e+00.

SNM/FB converges in 13 steps. At the final step, ||x¹³ - || » 3.0e-05, ||Y_FB(x¹³)|| » 6.2e-10, det L₁₃» -8.3e-05 (which indicates degeneracy). The rate of convergence is linear, with the ratio approaching 1/2.

The behavior of GNM/AS is reported in Table 2, and it clearly shows fast quadratic convergence.

The next four problems, taken from [14, Example 1-4] (Example 3.2 is slightly modified), are the KKT systems of the form (1.2), with g(z) = f¢(z), z Î ^p, where the objective function f : ^p ® will be specified for each example below.

Example 3.2. Let p = m = 2, f(z) = (z₁ + z₂)²/2 + (z₁ + z₂)³/3, G(z) = (z₁, z₂), z Î ², = 0, = 0. Semistability holds here, but for Y_NR (and hence, for Y_FB), BD-regularity is violated. The starting point is z⁰ = (1, 2), µ⁰ = (0.01, 0.01), with ||x⁰ - || » 2.2e+00, ||Y_FB(x⁰)|| » 1.7e+01, det L₀» 4.3e-04.

The behavior of SNM/FB is as follows: det L₁» 7.3e-10, det L₂» 0, but the corresponding linear system is solvable, and the method manages to escape the ''bad'' region. Specifically, det L₃» 4.3e-16, while det L₄» 2.5e-00, and the algorithm converges in 7 iterations. At the final step, ||x⁷ - || » 1.0e-16, ||Y_FB(x⁷)|| » 1.4e-16, det L₇» 2. The rate of convergence is superlinear.

The behavior of GNM/AS is reported in Table 3, and it also shows the superlinear rate.

Note that while SNM/FB and GNM/AS exhibit similar convergence for this problem, the performance of SNM/FB clearly depends on the specific implementation. In particular, solution of a given degenerate linear system depends on the linear solver, and can affect the overall convergence. Also, in the cases where BD-regularity is violated, different procedures to compute L_k could result (in general, not in this example) in different linear systems and some of them can be ill-conditioned close to the solution, preventing fast convergence of SNM/FB.

The next example shows that both 2-regularity of Y_S and weak regularity are important for fast convergence of GNM/AS.

Example 3.3. Let p = m = 2, f(z) = /2 + /3, G(z) = (z₁ - /2, z₁ + /2), z Î 2, = 0, = 0. Semistability is violated (and hence, BD-regularity for Y_NR and Y_FB is violated). For Y_S, 2-regularity holds, but weak regularity does not. The starting point is z⁰ = (0.1, 0.1), µ⁰ = (0.1, 0.1), with ||x⁰ - || » 2.0e-01, ||Y_FB(x⁰)|| » 1.3e-01, det L₀» 2.0e-01.

Both SNM/FB and GNM/AS converge in 12 steps, and ||x¹² - || » 2.4e-05, ||Y_FB(x¹²)|| » 8.4e-10, det L₁₂» 2.9e-04. The rate of convergence is linear with ratio 1/2.

Example 3.4. Let p = 2, m = 3, f(z) = z₁ + (+)/2, G(z) = (z₁, z₂, z₁ + z₂), z Î ², = 0, = (1, 0, 0). Semistability holds here, but for Y_NR (and hence, for Y_FB), BD-regularity is violated. The starting point is z⁰ = 0, µ⁰ = (1, 0.01, 0.01), with ||x⁰ - || » 2.2e-01, ||Y_FB(x⁰)|| » 1.2e-01, det L₀ = 0.

Though L₀ is singular, the SNM/FB linear system is solvable, and the method manages to escape. Specifically, det L₁ = -2, and the algorithm converges in 2 iterations: the second step gives the exact solution with det L₂» -3.5e-01. Note however that, as already discussed above, this behavior cannot be guaranteed for some other implementation of the linear system solver (and in general, of the procedure to compute L_k).

GNM/AS terminates after 1 step at the exact solution. The reason for this is that A₀ = {4, 5} , while , , and coincide with the corresponding components of . The iteration of GNM/AS terminates with the exact solution at the identification phase.

Example 3.5. Let p = m = 1, f(z) = z⁴/4, G(z) = z, z Î , = 0, = 0. Weak regularity holds here, but semistability (and hence, BD-regularity for Y_NR and Y_FB) and even 2-regularity for Y_S, are violated. The starting point is z⁰ = 1, µ⁰ = 0.1, with ||x⁰ - || » 1.0e+00, ||Y_FB(x⁰)|| » 9.1e-01, det L₀» 2.7e+00.

SNM/FB converges in 18 steps. At the final step, ||x¹⁸ - || » 6.8e-04, ||Y_FB(x¹⁸)|| » 3.1e-10, det L₁₃» 3.1e-06. The rate of convergence is linear with ratio approaching 2/3.

The behavior of GNM/AS is reported in Table 4, and it shows fast quadratic convergence.

The final example NCP taken from [15, Example 2.1].

Example 3.6. Let n = 2, F(x) = (-x₁ + x₂, -x₂), x Î ², = 0. BD-regularity holds for Y_NR (and hence, semistability also holds), but not for Y_FB. The starting point is x⁰ = (2, 4), with ||x⁰ - || » 4.5e+00, ||Y_FB(x⁰)|| » 5.8e+00, det L₀ = 0.

Here, SNM/FB fails to make a step. At the same time, GNM/AS terminates after 1 step at the exact solution. The reason for this is that A₀ = {1, 2} . Thus, the iteration of GNM/AS reduces to identifying the index sets.

Note that the problem in Example 3.6 is actually a linear complementarity problem, that is, NCP with affine F. We point out that in the case of affine F, just one step of GNM/AS gives the exact solution, provided the index sets are correctly identified. For example, this behavior is observed also for the problem badfree from the MCPLIB collection: once x^k is close to the solution = (0, 0, 0.5, 0.5, 1), GNM/AS produces x^k+1 = . At the same time, for x^k close to , a degenerate L_k is computed, and SNM/FB fails.

Acknowledgments. Research of the first two authors is supported by the Russian Foundation for Basic Research Grant 01-01-00810 and the Council for the State Support of Leading Scientific Schools Grant 00-15-96141. The third author is supported in part by CNPq Grants 300734/95-6(RN) and 471780/2003-0, by PRONEX-Optimization and by FAPERJ.

Received: 18/VIII/04. Accepted: 06/XII/04

#612/04.

• [1] S.C. Billups, Algorithms for complementarity problems and generalized equations. Technical Report 95-14, Computer Sciences Department, University of Wisconsin, Madison, Wisconsin, August 1995. PhD thesis.
• [2] J.F. Bonnans, Local analysis of Newton-type methods for variational inequalities and nonlinear programming, Applied Mathematics and Optimization, 29 (1994), 161-196.
• [3] A.N. Daryina, A.F. Izmailov and M.V. Solodov, A class of active-set Newton methods for mixed complementarity problems. IMPA preprint A256, October 2003 (revised April 2004). Accepted for publication in SIAM J. Optimization, 15 (2004), 409-429.
• [4] S.P. Dirkse and M.C. Ferris, MCPLIB: A collection of nonlinear mixing complementarity problems, Optimization Methods and Software, 5 (1995), 319-345.
• [5] A.S. El-Bakry, R.A. Tapia and Y. Zhang, A study of indicators for identifying zero variables in interior-point methods, SIAM Review, 36 (1994), 45-72.
• [6] F. Facchinei, A. Fischer and C. Kanzow, On the accurate identification of active non-straints, SIAM Journal on Optimization, 9 (1999), 14-32.
• [7] F. Facchinei, A. Fischer and C. Kanzow, On the identification of zero variables in an interior-point framework, SIAM Journal on Optimization, 10 (2000), 1058-1078.
• [8] F. Facchinei and S. Lucidi, Quadratically and superlinearly convergent algorithms for the solution of inequality constrained minimization problems, Journal of Optimization Theory and Applications, 85 (1995), 265-289.
• [9] F. Facchinei and J.-S. Pang, Finite-Dimensional Variational Inequalities and Complementarity Problems, Springer-Verlag, New York, NY, (2003).
• [10] M.C. Ferris, C. Kanzow and T.S. Munson, Feasible descent algorithms for mixed complementarity problems, Mathematical Programming, 86 (1999), 475-497.
• [11] M.C. Ferris and J.-S. Pang, Engineering and economic applications of complementarity problems, SIAM Review, 39 (1997), 669-713.
• [12] A.F. Izmailov and M.V. Solodov, Error bounds for 2-regular mappings with Lipschitzian derivatives and their applications, Mathematical Programming, 89 (2001), 413-435.
• [13] A.F. Izmailov and M.V. Solodov, The theory of 2-regularity for mappings with Lipschitzian derivatives and its applications to optimality conditions, Mathematics of Operations Research, 27 (2002), 614-635.
• [14] A.F. Izmailov and M.V. Solodov, Karush-Kuhn-Tucker systems: regularity conditions, error bounds and a class of Newton-type methods, Mathematical Programming, 95 (2003), 631-650.
• [15] T. De Luca, F. Facchinei and C. Kanzow, A theoretical and numerical comparison of some semismooth algorithms for complementarity problems, Computational Optimization and Applications, 16 (2000), 173-205.
• [16] T. Munson, F. Facchinei, M. Ferris, A. Fischer and C. Kanzow, The semismooth algorithm for large scale complementarity problems, INFORMS Journal of Computing, 13 (2001), 294-311.
• [17] S.J. Wright, Identifiable surfaces in constrained optimization, SIAM Journal on Control and Optimization, 31 (1993), 1063-1079.

Publication Dates

• Publication in this collection
  07 Nov 2005
• Date of issue
  Aug 2005

History

• Accepted
  06 Dec 2004
• Received
  18 Aug 2004