S Adly

A L Dontchev

M Théra

ON ONE-SIDED LIPSCHITZ STABILITY OF SET-VALUED CONTRACTIONS

Keywords: Composition, Differential inclusions, Fixed points, Lyusternik-Graves theorem, Lipschitz stability, Set-valued mappings

We give conditions under which the distance from a point x to the set of fixed points of the composition of the set-valued mappings F and G is bounded by a constant times the smallest distance between F -1 (x) and G (x). This estimate allows us to significantly sharpen a result by T.-C. Lim [10] regarding fixed-points stability of set-valued contractions. A global version of the Lyusternik-Graves theorem is obtained from this estimate as well. We apply the generalization of Lim's result to establish one-sided Lipschitz properties of the solution mapping of a differential inclusion with a parameter.

In a path-breaking paper, Teck-Cheong Lim [START_REF] Lim | On fixed-point stability for set-valued contractive mappings with applications to generalized differential equations[END_REF] showed that the Pompeiu-Hausdorff distance between the sets of fixed points of two setvalued contraction mappings is bounded by a constant times the uniform Pompeiu-Hausdorff distance between the mappings. Specifically, he proved the following result:

Lim's Lemma (Lim [10, Lemma 1]). Let X be a complete metric space and let T 1 and T 2 map X into the family of nonempty closed subsets of X . Suppose that both T 1 and T 2 are Lipschitz continuous on X with the same Lipschitz constant 1 ∈ [0, 1). Then

haus (Fix(T 1 ), Fix(T 2 )) ≤ 1 1 - sup x∈X haus (T 1 (x), T 2 (x)) (1)
Before going further, let us fix the notation and terminology. Throughout this article, X and Y are metric spaces unless stated otherwise.

(D) = x ∈ X | ∃ y ∈ Y such that (x, y) ∈ D .
Any mapping, possibly set-valued and acting from a space X to (the subsets of) a space Y , is denoted by

F : X ⇒ Y . The domain of F is dom F = x ∈ X | F (x) = ∅ , its range is rge F = y | ∃ x with y ∈ F (x) , and its graph is gph F = (x, y) ∈ X × Y | y ∈ F (x) . The inverse F -1 : Y ⇒ X of a mapping F : X ⇒ Y is defined as F -1 (y) = x | y ∈ F (x)
. By a composition of a mappings B : Z ⇒ Y with a mapping A : X × P ⇒ Z we mean a mapping C = B • A : X × P ⇒ Y whose graph contains all points (x, p, y) ∈ gph C for which there exists z ∈ Z such that (x, p, z) ∈ gph A and (z, y) ∈ gph B; in other words, C (x, p) = y ∈ B(z) | z ∈ A(x, p) . We denote the set of fixed points of a mapping F : X ⇒ X by Fix(F ); that is, Fix(F

) := x ∈ X | x ∈ F (x) Note that for mappings F : Y ⇒ X and G : X ⇒ Y we have Fix(F • G ) = x ∈ X : F -1 (x) ∩ G (x) = ∅ . A set-valued mapping F : X ⇒ Y is said to be Lipschitz continuous relative to a set D ⊂ dom F when F is closed valued 1 on D and there exists a constant ≥ 0 (Lipschitz constant) such that haus (F (x), F (x )) ≤ (x, x ) for all x, x ∈ D (2)
When D = X we say that F is Lipschitz continuous. In that case, dom F = X . If a mapping F acting from a metric space X to a complete metric space Y is Lipschitz continuous, then F necessarily has a closed graph. Indeed, let (x k , y k ) ∈ gph F and (x k , y k ) → (x, y); then

d(y, F (x)) ≤ (y, y k ) + d(y k , F (x k )) + haus (F (x k ), F (x)) ≤ (y, y k ) + (x k , x) → 0,
and, hence, (x, y) ∈ gph F . It is well documented in the literature that the Lipschitz continuity of a mapping F is equivalent to two properties of its inverse F -1 : metric regularity and linear openness. Without going into this further, we use the fact established in [START_REF] Dontchev | Implicit Functions and Solution Mappings[END_REF]Proposition 3C.1] that a mapping G : Y ⇒ X is Lipschitz continuous with constant if and only if

d(x, G (y)) ≤ d(y, G -1 (x)) for every x ∈ X , y ∈ Y (3)
Going back to Lim's lemma [START_REF] Lim | On fixed-point stability for set-valued contractive mappings with applications to generalized differential equations[END_REF] we started with, first observe that it implies the standard Nadler fixed point theorem for set-valued mappings [START_REF] Nadler | Multi-valued contraction mappings[END_REF]; indeed, if T 1 = T 2 = T then the right side of (1) is zero and, hence, so is the left side, thus Fix(T ) is nonempty. Right after the statement of his lemma, Lim stated a theorem (Theorem 1) which follows immediately from the lemma and which Lim apparently considered as the central result in his paper. This theorem says that when T n is a sequence of set-valued mappings T n : X ⇒ X each of which is Lipschitz continuous with the same Lipschitz constant ∈ [0, 1), then

lim n→∞ sup x∈X haus (T n (x), T 0 (x)) = 0 ⇒ lim n→∞ haus (Fix(T n ), Fix(T 0 )) = 0 (4)
In Theorems 2 and 3 in [START_REF] Lim | On fixed-point stability for set-valued contractive mappings with applications to generalized differential equations[END_REF], Lim presented extensions of this implication in some particular spaces and compact-valued nonexpansive mappings. Finally, he applied his Theorem 1 to show continuity of the solution mapping for a differential inclusion with respect to the initial conditions.

Observe that Lim's lemma is a quantitative result and, therefore, it is a natural tool to establish quantitative stability of perturbed fixed points. Specifically, the following corollary of Lim's lemma shows Lipschitz continuity of fixed points: Theorem 2 (Lipschitz continuity of fixed points). Let X be a complete metric space and P be a metric space. Consider a mapping M : P × X ⇒ X with the following properties:

(i) M (p, •) is Lipschitz continuous with a Lipschitz constant ∈ [0, 1) uniformly in p ∈ P ; (ii) M (•, x) is Lipschitz continuous with a Lipschitz constant uniformly in x ∈ X .
Then the mapping p → Fix(M (p, •)) is Lipschitz continuous with a constant /(1 -).

Proof. We will first show that for each p ∈ P the set Fix(M (p,

•)) is closed. Let u k ∈ Fix(M (p, •)) for all k and u k → u. Then u k ∈ M (p, u k ) for all k. We have d(u, M (p, u)) ≤ d(u, M (p, u k )) + haus (M (p, u k ), M (p, u)) ≤ (u, u k ) + (u, u k ) → 0 as k → ∞ Hence, d(u, M (p, u)) = 0, that is, u ∈ M (p, u), which implies that u ∈ Fix(M (p, •)), thus Fix(M (p, •)) is closed. Let p 1 , p 2 ∈ P . Applying Lim's lemma with T 1 (x) = M (p 1 , x) and T 2 (x) = M (p 2 , x) we have haus Fix(M (p 1 , •)), Fix(M (p 2 , •)) ≤ 1 1 - sup x∈X haus (M (p 1 , x), M (p 2 , x)) ≤ 1 - (p 1 , p 2 )
As a corollary of our main result presented in Theorem 6, we will obtain the following generalization of Lim's lemma, where the Pompeiu-Hausdorff distance is replaced by either of the excesses.

Theorem 3 (generalized Lim's lemma). On the conditions of Lim's lemma, e(Fix(T

i ), Fix(T j )) ≤ 1 1 - sup x∈X e(T i (x), T j (x)) fori, j = 1, 2
Theorem 3 allows us to obtain results that are parallel to Theorem 2 but for the weaker one-sided Lipschitz continuity property. First, we need the following definition:

Definition. A mapping F : X ⇒ Y is said to be outer [resp. inner] Lipschitz continuous at x relative to a set D ⊂ dom F when x ∈ D and there exists a constant ≥ 0 such that e(F (x), F ( x)) ≤ (x, x) [resp. e(F ( x), F (x)) ≤ (x, x)] for all x ∈ D (5)
When D = X they we say that F is outer [resp. inner] Lipschitz continuous at x. In that case, dom F = X .

Observe that both the outer and the inner Lipschitz continuity are much weaker properties than the Lipschitz continuity. Directly from the definition, when F is both outer and inner Lipschitz continuous at a point, then it is Lipschitz continuous. The property of outer Lipschitz continuity was introduced by Robinson [START_REF] Robinson | Generalized equations and their solutions. I. Basic theory. Point-to-set maps and mathematical programming[END_REF] under the name "upper Lipschitz continuity," the main motivation being the important result due to him that every mapping F : n ⇒ m whose graph is the union of finitely many polyhedral convex sets is outer Lipschitz continuous at any point of its domain.

From Theorem 3 we can obtain the following result, which strengthens Theorem 2.

Theorem 4. Let X be a complete metric space and P be a metric space. Consider a mapping M : P × X ⇒ X with the following properties:

(i) M (p, •) is Lipschitz continuous with a Lipschitz constant ∈ [0, 1) uniformly in p ∈ P ; (ii) M (•, x) is outer [resp. inner] Lipschitz continuous at p with a constant uniformly in x ∈ X . Then the mapping p → Fix(M (p, •)) is outer [resp. inner] Lipschitz continuous at p with a constant /(1 -). Proof. Applying Theorem 3 with T 1 (x) = M (p, x) and T 2 (x) = M ( p, x) we have e(Fix(M (p, •), Fix(M ( p, •)) ≤ 1 1 - sup x∈X e(M (p, x), M ( p, x)) ≤ 1 - (p, p)
Arutyunov derived in [START_REF] Arutyunov | Covering mappings in metric spaces and fixed points[END_REF] a coincidence theorem which we present here in the following form, where the assumptions about completeness/closedness of graphs are different from those originally used: Theorem 2]). Consider mappings F : Y ⇒ X and G : X ⇒ Y and assume that one of the set gph F , gph G is complete and the other is closed. Let and be nonnegative constants such that < 1. Suppose that F is Lipschitz continuous with constant and G is Lipschitz continuous with constant . Then for every > 0 and every x ∈ X there exists a point ∈ X such that

Theorem 5 (Arutyunov [1,
F -1 ( ) ∩ G ( ) = ∅ and (x, ) ≤ 1 - d(F -1 (x), G (x)) +
The main part of the result given next as Theorem 6; namely, the implication (i) ⇒ (ii), is equivalent, up to some adjustments regarding closedness/completeness of sets, to Arutyunov's theorem.

Theorem 6 (estimate for fixed points of composition). Let and be nonnegative constants such that < 1. Consider a mapping F : Y ⇒ X with closed and nonempty graph and P X (gph F ) complete. Then the following are equivalent:

(i) The mapping F is Lipschitz continuous with constant ;

(ii) dom F = Y and for any mapping G : X ⇒ Y with closed graph which is Lipschitz continuous with constant and such that P Y (gph G ) is complete, the following inequality holds:

d(x, Fix(F • G )) ≤ 1 - d(F -1 (x), G (x)) for everyx ∈ X (6)
Proof. Let (i) hold, let G be as required in (ii) and let x ∈ X . If x rge F or G (x) = ∅, the right side of ( 6) is +∞ and we are done. Otherwise, let y ∈ G (x). Choose > 0. Since by definition dom

F = Y , there exists u ∈ F (y) such that (x, u) < d(x, F (y)) + If u = x then x ∈ Fix(F • G )
and the left side of ( 6) is zero, so there is nothing more to prove. Otherwise, we have

d(x, F (y)) < (1 + ) (u, x) (7)
Denote a := (1 + ) (u, x); then

a ≤ (1 + )(d(x, F (y)) + ) (8)
We will construct a sequence (x k , y k ) with the following properties:

y k ∈ F -1 (x k+1 ) ∩ G (x k ) (9)
and

(y k , y k-1 ) < a( ) k-1 , (x k+1 , x k ) < a( ) k ( 10
)
Let x 0 = x and y 0 = y. In the first step, observe that, from [START_REF] Dontchev | Implicit Functions and Solution Mappings[END_REF],

d(x 0 , F (y 0 )) < a
Then there exists x 1 ∈ F (y 0 ) with

(x 1 , x 0 ) < a( ) 0
Furthermore, since d(y 0 , G (x 1 )) ≤ haus (G (x 0 ), G (x 1 )) ≤ (x 0 , x 1 ) < a, there exists y 1 ∈ G (x 1 ) such that (y 0 , y 1 ) < a( ) 0

From the Lipschitz continuity of F , we have d(x 1 , F (y 1 )) ≤ haus ((y 0 ), F (y 1 )) ≤ (y 0 , y 1 ) < a( ) 0 = a( ) Hence, there exists x 2 ∈ F (y 1 ) with

(x 2 , x 1 ) < a( ) 1
We obtain ( 9) and ( 10) for k = 1.

Proceeding by induction, suppose that we have already found a sequence (x k , y k ) satisfying ( 9) and ( 10) for k = 1, , j , for some j > 1. Then

d(y j , G (x j +1 )) ≤ haus (G (x j ), G (x j +1 )) ≤ (x j , x j +1 ) < a( ) j
Hence, one can find y j +1 ∈ G (x j +1 ) such that (y j , y j +1 ) < a( ) j

The Lipschitz continuity of F gives us d(x j +1 , F (y j +1 )) ≤ haus (F (y j ), F (y j +1 )) ≤ (y j , y j +1 ) < a( ) j = a( ) j +1 Then there exists x j +2 ∈ F (y j +1 ) with

(x j +2 , x j +1 ) < a( ) j +1
The induction step is complete.

For k > m + 1 > 1, from [START_REF] Lim | On fixed-point stability for set-valued contractive mappings with applications to generalized differential equations[END_REF] we have

(x k , x m ) ≤ k-1 i=m (x i+1 , x i ) < a( ) m 1 - and (y k , y m ) ≤ k-1 i=m (y i+1 , y i ) < a( ) m 1 -
Hence, (x k , y k ) is a Cauchy sequence. Since (x k , y k ) belongs to the complete set P X (gph F ) × P Y (gph G ) for all k we conclude that this sequence is convergent to some ( x, ŷ). Since both F and G have closed graphs, we have x ∈ F (ŷ) and ŷ ∈ G ( x), hence x ∈ Fix(F • G ).

Note that

(x k , x) ≤ k-1 i=0 (x i+1 , x i ) < a k-1 i=0 ( ) i < a (1 -)
Passing to the limit with k → ∞, from (8) we obtain

d(x, Fix(F • G )) ≤ ( x, x) ≤ (1 + )(d(x, F (y)) + ) 1 -
Since can be arbitrarily small, we get

d(x, Fix(F • G )) ≤ 1 1 - d(x, F (y)) (11)
This combined with (3) gives us

d(x, Fix(F • G )) ≤ 1 - d(y, F -1 (x))
Taking into account that y can be arbitrarily chosen in G (x) we obtain [START_REF] Cernea | Stability of solution sets of nonlinear integrodifferential inclusions[END_REF]. Now, suppose that (ii) is satisfied. Let y, y ∈ Y and let x ∈ F (y). Choose G (x) = y for all x ∈ X . Clearly, P Y (gph G ) is complete and G is Lipschitz continuous with constant = 0; moreover, d(x, Fix(F • G )) = d(x, F (y )). From (6),

d(x, F (y )) ≤ d(y , F -1 (x)) ≤ (y, y )
Taking supremum with respect to x ∈ F (y) in the left side of this inequality, we obtain that F is Lipschitz continuous.

Remark 1. Note that the estimate ( 6) is sharp in the sense that if the left side is zero, so is the right side.

Remark 2. One of the referees of the first version of this paper suggested the following proof that (a part of) Theorem 6 follows from Theorem 5. Put U = P X (gph F ) and apply Theorem 5 with F being the identity mapping from U to U and G identified with the composition F • G restricted on U , which is Lipschitz continuous with Lipschitz constant < 1. For any x ∈ U from Theorem 5 we have

d(x, Fix(F • G )) ≤ 1 1 - d(x, (F • G )(x))
Further, for any y ∈ G (x) we have

d(x, Fix(F • G )) ≤ d(x, F (y))
Combining the last two estimates with (3) and noting that y is an arbitrary point in G (x), we obtain [START_REF] Azé | On the dependence of fixed point sets of pseudo-contractive multifunctions. Application to differential inclusions[END_REF]. When x U the right side of ( 5) is +∞ and there is nothing more to prove.

Remark 3. Lim proved his lemma in [START_REF] Lim | On fixed-point stability for set-valued contractive mappings with applications to generalized differential equations[END_REF] by using an iterative procedure which is quite different from the one used in the proof of Theorem 6.

In particular, this procedure involves a construction of just one sequence of iterates x k for both mappings T 1 and T 2 . Also, the iteration we use is related but different from the one in [START_REF] Arutyunov | Covering mappings in metric spaces and fixed points[END_REF].

Proof of Theorem 3. Apply Theorem 6 with X = Y , F = I , the identity mapping, and G = T 1 . Then G is Lipschitz continuous with Lipschitz constant ∈ [0, 1) and from [START_REF] Cernea | Stability of solution sets of nonlinear integrodifferential inclusions[END_REF] we have that for any x ∈ X ,

d(x, Fix(T 1 )) ≤ 1 1 - d(x, T 1 (x)) (12)
Let T 2 : X ⇒ Y be Lipschitz continuous with Lipschitz constant . Clearly, (12) holds for every x ∈ Fix(T 2 ) and, hence, taking supremum with respect to such x, we get

e(Fix(T 2 ), Fix(T 1 )) ≤ 1 1 - sup x∈T 2 (x) d(x, T 1 (x)) ≤ sup x∈X 1 1 - e(T 2 (x), T 1 (x))
By symmetry, the proof is complete.

Proof of Theorem 5. It is easy to modify the proof of Theorem 6 for the case when one of the sets gph F , gph G is complete while the other is closed. Indeed, if gph F is complete, then (x k , y k ) is convergent, and hence the limit ( x, ŷ) satisfies x ∈ F (ŷ) and ŷ ∈ G ( x), Let x ∈ X and > 0. From (6), there exists

∈ (F • G )( ), or equivalently, F -1 ( ) ∩ G ( ) = ∅, such that (x, ) ≤ 1 - d(F -1 (x), G (x)) +
Another application of Theorem 6 is the following result which is a version of a general paradigm in nonlinear analysis the most popular form of which appears in the basic inverse function theorem, for more see [7, chap. 5]. A local version of it stems from early works of Lyusternik and Graves and is most known as the Lyusternik-Graves theorem. Here we Let a > 0 and let B be a bounded set in n . Consider a mapping R acting from B into the family of nonempty closed subsets of B which is Lipschitz continuous with respect to the Pompeiu-Hausdorff metric with a Lipschitz constant ; that is, haus (R (x), R (y)) ≤ xy for all x, y ∈

Given an initial point b ∈ B, consider the initial-value problem

ẋ(t ) ∈ R (x(t )) for a.e. t ∈ [0, a], x(0) = b ( 14
)
By a solution of ( 14), we mean any absolutely continuous function which satisfies [START_REF] Yin | Metric regularity of composite multifunctions in Banach spaces[END_REF] for almost every (a.e.) t ∈ [0, a], with x( 0

such that y(t ) = b + t 0 z( )d for all t ∈ [0, a] ,
where the integral is in the sense of Lebesgue. We will first show that T (b, x) = ∅ when x ∈ C ([0, a], B) Choose any x ∈ C ([0, a], B) For each t ∈ [0, a], let z(t ) be a projection of 0 ∈ n on the closed set R (x(t )). Since both x and R are continuous, we obtain that the function z defined in such a way is measurable; since it is uniformly bounded, it is integrable. Then y ∈ T (b, x) when y(t

) = b + t 0 z( )d for t ∈ [0, a].
To prove Lipschitz continuity of T (b, •) we use a similar argument. Let x 1 , x 2 ∈ C ([0, a], B) and let y 1 ∈ T (b, x 1 ). From the very definition of the mapping T , there exists a measurable selection z 1 of R (x 1 ) such that y 1 (t ) = b + t 0 z 1 ( )d for all t ∈ [0, a]. For each t ∈ [0, a], let z 2 (t ) be a projection of z 1 (t ) on R (x 2 (t )). Since z 1 is measurable and R is continuous, we obtain that the function z 2 defined in such a way is measurable; since it is uniformly bounded, it is integrable. Let y 2 (t ) = b + t 0 z 2 ( )d for t ∈ [0, a]. Then, from the Lipschitz continuity property (13) of R we obtain

y 1 -y 2 C ≤ a 0 z 1 ( ) -z 2 ( ) d ≤ a 0 x 1 ( ) -x 2 ( ) d ≤ a x 1 -x 2 C Hence, the mapping T (b, •) is Lipschitz continuous with a Lipschitz constant a which is independent of b. Let b 1 , b 2 ∈ B and denote T 1 (x) = T (b 1 , x) and T 2 (x) = T (b 2 , x). Observe that haus (T 1 (x), T 2 (x)) ≤ b 1 -b 2 (15)
Indeed, let y 1 ∈ T 1 (x). Then there exists an associate selection z of R (x)

such y 1 (t ) = b 1 + t 0 z( )d for all t ∈ [0, a]. Clearly y 2 defined as y 2 (t ) = b 2 + t 0 z( )d , t ∈ [0, a], is from T 2 (x)
. This yields (15). Denote by S (b) the set of solutions of the differential inclusion [START_REF] Yin | Metric regularity of composite multifunctions in Banach spaces[END_REF] for a fixed initial condition b ∈ B. Then S (b) = Fix(T (b, •)). From Lim's lemma we obtain the following theorem: Theorem 8. Consider the differential inclusion [START_REF] Yin | Metric regularity of composite multifunctions in Banach spaces[END_REF] for a Lipschitz continuous mapping R with a Lipschitz constant as in [START_REF] Singh | A note on fixed point stability for generalized set-valued contractions[END_REF]. Choose a > 0 such that a < 1. Then the solution mapping S of ( 14) is Lipschitz continuous from the set B to the collection of nonempty subsets of functions from the space C ([0, a], B). In Lim's paper [START_REF] Lim | On fixed-point stability for set-valued contractive mappings with applications to generalized differential equations[END_REF], it is assumed that R is in addition convex-valued and the conclusion is that the solution mapping S is merely continuous, obtained as an application of [10, Theorem 1].

We can obtain even more general results for a differential inclusion with a parameter, that is,

ẋ(t ) ∈ R (p, x(t )) for a.e. t ∈ [0, a], x(0) = b, ( 16
)
where now b ∈ B is fixed and p is a parameter from a metric space P with a reference value p. By repeating the above argument where we use Theorem 4 instead of Lim's lemma, we obtain: Theorem 9 can be easily extended to a nonautonomous setting, as in Lim's original paper [START_REF] Lim | On fixed-point stability for set-valued contractive mappings with applications to generalized differential equations[END_REF], or even in a much more general context of infinite-dimensional spaces. Specifically, in a series of papers A. Cernea considered various nonlinear differential inclusions and proved stability of their solutions by using Lim's lemma. In particular, in [START_REF] Cernea | Stability of solution sets of nonlinear integrodifferential inclusions[END_REF] he considered nonlinear integrodifferential inclusions of the following form: ẋ(t ) ∈ F (t , x(t ), V (x)(t )) for a.e. t ∈ [0, ∞), x(0) = x 0 , (

where F : [0, +∞) × X × X ⇒ X , X is a separable Banach space and V : C ([0, +∞), X ) → C ([0, +∞), X ) is the nonlinear Volterra integral operator

V (x)(t ) = t 0 k(s, x(s))ds,
where k : [0, +∞) × X → X is Lipschitz continuous. Applying Lim's lemma, Cernea proved Lipschitz continuity of the Lipschitz mapping on an infinite time horizon. To do that, instead of considering sufficiently small time interval as in Theorem 8 above, he used a weighted L 2 space with a sufficiently large weight to obtain Lipschitz continuity of the mapping defining the inclusion (17) with a Lipschitz constant less than 1. Then he followed the approach based on Lim's lemma. By using the extension of Lim's lemma in this article, this result can be extended to differential inclusion of the form (17) involving a parameter p ∈ P and employing, as in Theorem 9, the outer [resp. inner] Lipschitz continuity.

  Any metric is denoted by (•, •), any norm by • , and d(x, C ) is the distance from a point x ∈ X to a subset C ⊂ X ; that is d(x, C ) = inf (x, x ) x ∈ C whenever C = ∅ and d(x, ∅) = ∞. The closed ball centered at x with radius r is denoted by r (x), and the closed unit ball is . The excess from a set C to a set D is e(C , D) = sup x∈C d(x, D) under the convention e(∅, D) = 0 for D = ∅ and e(D, ∅) = +∞ for any D; thus, if e(C , D) < ∞ then D must be nonempty. The Pompeiu-Hausdorff distance between C and D is haus (C , D) = max e(C , D), e(D, C ) . The smallest distance between two sets C and D is denoted by d(C , D), that is, d(C , D) = inf (x , x) x ∈ C , x ∈ D . If one of the sets C and D is empty, we set d(C , D) = haus (C , D) = +∞. Given a set D ⊂ X × Y we denote by P X (D) the projection of D on X , that is, P X

  ) = b. Let C ([0, a], B) be the space of continuous functions on the interval [0, a] with values in B, equipped with the usual supremum norm • C . Define the mapping C ([0, a], B) x → T (b, x) = y ∈ C ([0, a], B) there exists a Lebesgue measurable function z such that z( ) ∈ R (x( )) a.e. ∈ [0, a]

Theorem 9 .

 9 Consider the differential inclusion (16), where the mapping R is Lipschitz continuous with respect to x ∈ B with constant , uniformly in p ∈ P , and outer [resp. inner] Lipschitz continuous with respect to p at p with constant , uniformly in x ∈ B. Let a < 1. Then the solution mapping P p → S (p) = x ∈ C ([0, a], B) | x solves (16) for p ∈ P is outer [resp. inner] Lipschitz continuous at p with a constant /(1a).

Following[START_REF] Dontchev | Implicit Functions and Solution Mappings[END_REF], we require that Lipschitz continuous mappings must be closed valued to avoid the situation when Lipschitz continuity does not entail continuity. This is also assumed implicitly in[START_REF] Lim | On fixed-point stability for set-valued contractive mappings with applications to generalized differential equations[END_REF].

ACKNOWLEDGMENTS

The authors wish to thank the anonymous referees for their valuable comments and suggestions.

FUNDING A. L. Dontchev was supported by the National Science Foundation, Grant DMS 1008341 through the University of Michigan. M. Théra was Supported by Grant DP 11010211 of the Australian Research Council through the University of Ballard.

deduce from Theorem 6 the following global version of the Lyusternik-Graves theorem:

Theorem 7 (Global Lyusternik-Graves). Let X be a complete metric space and Y be a Banach space. Consider mappings : X ⇒ Y and : X ⇒ Y and suppose that both -1 and are globally Lipschitz continuous with constant > 0 and > 0, respectively. If

Applying Theorem 6 with F = -1 and G (•) = -(•) + y to the last inequality, from (6) we get

There is a vast literature on fixed points of compositions which, however, cover topics that are different from the subject of the present paper. Results related to the ones obtained in the current paper are given in [2-5, 8, 9, 13, 14].

In his paper, Lim [START_REF] Lim | On fixed-point stability for set-valued contractive mappings with applications to generalized differential equations[END_REF]Theorem 4] presented an application of his Theorem 1 given by ( 4) to obtain a stability result for the solution mapping of an initial value problem for a differential inclusion. Specifically, he proved continuity of the set of solutions, in the space of the continuous functions equipped with the supremum norm, with respect to the initial condition. However, Lim did not mention anything about Lipschitz continuity of the solution set although such a statement easily follows from his analysis. Also, some of his assumptions he maid are not necessary to obtain such a result. In further lines we will first present the setting of Lim's result (Theorem 4) with some adjustments in terminology and notation, and then show how to obtain Lipschitz continuity. Then we generalize this result to one-sided Lipschitz stability of the solution mapping with respect to a parameter. In order to simplify the presentation, we work in the Euclidean space n and consider an autonomous differential inclusion over a finite time interval.


Samir Adly
email: samir.adly@unilim.fr

Hadia Rammal
email: hadia.rammal@xlim.fr

A new method for solving Pareto eigenvalue complementarity problems

Keywords: Complementarity problems, One-constrained eigenvalue problems, Complementarity functions, Semismooth Newton Method, Lattice Projection Method, Bi-eigenvalue complementarity problems

In this paper, we introduce a new method, called the Lattice Projection Method (LPM), for solving eigenvalue complementarity problems. The original problem is reformulated to find the roots of a nonsmooth function. A semismooth Newton type method is then applied to approximate the eigenvalues and eigenvectors of the complementarity problems. The LPM is compared to SNM min and SNM FB , two methods widely discussed in the literature for solving nonlinear complementarity problems, by using the performance profiles as a comparing tool (Dolan, Moré in Math. Program. 91:201-213, 2002). The performance measures, used to analyze the three solvers on a set of matrices mostly taken from the Matrix Market (Boisvert et al. in The quality of numerical software: assessment and enhancement, pp. [125][126][127][128][129][130][131][132][133][134][135][136][137] 1997), are computing time, number of iterations, number of failures and maximum number of solutions found by each solver. The numerical experiments highlight the efficiency of the LPM and show that it is a promising method for solving eigenvalue complementarity problems. Finally, Pareto bi-eigenvalue complementarity problems were solved numerically as an application to confirm the efficiency of our method.

Introduction

Since the last three decades, the subject of complementarity problems has become one of the most well-established disciplines within mathematical programming. Due to its various applications in engineering, economics and sciences, it is not surprising that many researchers have focused their attention on developing both theoretical results and efficient numerical methods for solving these problems (see the exhaustive survey [START_REF] Ferris | Engineering and economic applications of complementarity problems[END_REF]). Among these problems, we focus, in our study, on the Eigenvalue Complementarity Problems (EiCP), also called, cone-constrained eigenvalue problems [1, 18-20, 32, 33, 36, 38, 39] that are considered as nonlinear optimization problems. These problems consist to find a scalar (eigenvalue) and a nonzero vector (eigenvector) satisfying a complementarity relation over a closed convex cone in a finite-dimensional space. This case recovers the classical eigenvalue problem when the closed convex cone coincides with the whole space. The Pareto-eigenvalue problems correspond to the case when the closed convex cone coincides with the positive orthant R n + . A wide variety of applications in sciences and engineering require the numerical computation of the solutions of EiCP such as the dynamic analysis of structural mechanical systems, vibro-acoustic systems, electrical circuit simulation, signal processing, fluid dynamic, contact problem in mechanics (see for instance [START_REF] Martins | Dynamic stability of finite dimensional linearly elastic systems with unilateral contact and Coulomb friction[END_REF][START_REF] Martins | Stability of finite-dimensional nonlinear elastic systems with unilateral contact and friction[END_REF][START_REF] Martins | Bifurcations and instabilities in frictional contact problems: theoretical relations, computational methods and numerical results[END_REF][START_REF] Martins | Computation of bifurcations and instabilities in some frictional contact problems[END_REF][START_REF] Martins | The directional instability problem in systems with frictional contacts[END_REF]). In some applications, the knowledge of eigenvalues can avoid the instability and unwanted resonance for a given system. In mechanical structures for example, eigenvalues are closely related to resonance frequency and to the stability analysis of the corresponding dynamical systems. The computation of eigenvalues becomes crucial in order to be located and damped eigenvalues corresponding to unstable modes or involving large vibrations (see [START_REF] Tisseur | The quadratic eigenvalue problem[END_REF]).

The theoretical spectral analysis is now well-developed (see [START_REF] Pinto Da Costa | Cone-constrained eigenvalue problems: theory and algorithms[END_REF]), while the construction of robust and efficient algorithms, for solving EiCP, needs further investigation. Note that the number of Pareto-eigenvalues grows exponentially with the dimension of the problem, for example, a matrix of order 25 may have more than 3 million Pareto-eigenvalues. Therefore, finding the whole Pareto-spectrum of a medium or a large size matrix is not an easy task. Recently, S. Adly and A. Seeger [START_REF] Adly | A nonsmooth algorithm for cone-constrained eigenvalue problems[END_REF] reformulated the EiCP as an equivalent nonsmooth system of equations, by using the Nonlinear Complementarity Functions (NCP-functions) φ min and φ FB , and then applied the corresponding Semismooth Newton Method (SNM min and SNM FB ). The major contribution of this paper is to introduce a new method, namely the Lattice Projection Method (LPM for short), for solving EiCP. This method is not based on the complementarity approach since it does not use any NCP-functions. The LPM is compared to SNM min and SNM FB by using the performance profiles as a comparison tool [START_REF] Dolan | Benchmarking optimization software with performance profiles[END_REF] and the numerical experiments highlight that the LPM solver is efficient and robust. We have also studied under which conditions all the elements of the Clarke generalized Jacobians ∂Φ(z * ) defined in [START_REF] Judice | The eigenvalue complementarity problem[END_REF] and ∂Φ LPM (z * ) defined in [START_REF] Martins | Stability of finite-dimensional nonlinear elastic systems with unilateral contact and friction[END_REF] at solutions z * and z * of ( 17) and [START_REF] Martins | Dynamic stability of finite dimensional linearly elastic systems with unilateral contact and Coulomb friction[END_REF], respectively, are nonsingular.

The structure of this paper is as follows. In the next section, some theoretical results related to EiCP are briefly discussed. In Sect. 3, we establish a short description of the Semismooth Newton Method by presenting the algorithm and its convergence theorem. Moreover, we reformulate the EiCP as an equivalent nonsmooth system of equations, by using the Nonlinear Complementarity Functions (NCP-functions) φ min and φ FB . The new method, called Lattice Projection Method (LPM), for solving EiCP is discussed in Sect. 4. In Sect. 5, we study the nonsingularity conditions of the Clarke generalized Jacobians ∂Φ(z * ) and ∂Φ LPM (z * ) defined in [START_REF] Judice | The eigenvalue complementarity problem[END_REF] and [START_REF] Martins | Stability of finite-dimensional nonlinear elastic systems with unilateral contact and friction[END_REF], respectively. Some examples are discussed to illustrate this regularity concept. We provide, in Sect. 6, some numerical results related to some matrices of order 3, 4 and 5 (known to have the maximum number of Pareto-eigenvalues). In Sect. 7, we give some extensive numerical tests showing that the LPM is extremely robust and very effective. Performance profiles and data, mostly taken from the Matrix Market [START_REF] Boisvert | Matrix market: a web resource for test matrix collections[END_REF], are used to analyze the performance of the three algorithms: LPM, SNM min and SNM FB . Computing time, number of iterations, number of failures and maximum number of eigenvalues found by each solver are used as performance measures to compare these algorithms. In Sect. 8, we show that the LPM can be generalized to solve the bivariate Eigenvalue Complementarity problems (bi-EiCP) [START_REF] Seeger | Eigenvalue analysis of equilibrium processes defined by linear complementarity conditions[END_REF]. Numerical results are then reported to confirm the efficiency and the robustness of the LPM. Finally, in Sect. 9, we close up the paper with a conclusion.

Pareto eigenvalue complementarity problems

For a real matrix A ∈ M n (R) of order n, the eigenvalue complementarity problem EiCP associated to the Pareto cone R n + is defined as follows:

find λ ∈ R and x ∈ R n \ {0}, such that x ≥ 0, λx -Ax ≥ 0, x, λx -Ax = 0. (1)
The scalar λ and the vector x satisfying system (1) are respectively called a Paretoeigenvalue of A and an associated Pareto-eigenvector. The set of all eigenvalues is called the Pareto-spectrum of A, and it is defined by

σ (A) = λ ∈ R : ∃x ∈ R n \ {0}, 0 ≤ x ⊥ (λx -Ax) ≥ 0 .
The following result, due to A. Seeger in [START_REF] Seeger | Eigenvalue analysis of equilibrium processes defined by linear complementarity conditions[END_REF], fully characterized the Pareto-spectrum of (1).

Lemma 1 [START_REF] Seeger | Eigenvalue analysis of equilibrium processes defined by linear complementarity conditions[END_REF] The scalar λ is a Pareto-eigenvalue of A if and only if there exist a nonempty index set I ⊆ N = {1, . . . , n}, and a vector η ∈ R |I | such that

A I η = λη, η ∈ int R |I | + , ( 2
)
j ∈I a ij η j ≤ 0, ∀i ∈ N \ I, ( 3
)
where A I ∈ M |I | (R) denotes the principal submatrix of A which is obtained by deleting the ith row and the ith column of A, whenever i / ∈ I . In this case, the vector x ∈ R n given by

x i = η i if i ∈ I, 0 if i ∈ N \ I.
is a Pareto-eigenvector of A associated to the eigenvalue λ.

Remark 1 This lemma shows that solving the EiCP is equivalent to solve 2 n -1 classical eigenvalue subproblems involving the sub-matrices A I and select only the solutions satisfying the condition (3). This result can be used to find all Pareto-eigenvalues of a matrix A of small order n.

The Pareto capacity of the space M n (R), is defined by

π n = max A∈M n (R) card σ (A) ,
where card[σ (A)] is the cardinality of the Pareto-spectrum of A. As proved in [START_REF] Pinto Da Costa | Cone-constrained eigenvalue problems: theory and algorithms[END_REF][START_REF] Seeger | On cardinality of Pareto spectra[END_REF], one has

3 2 n-1 -1 ≤ π n ≤ n2 n-1 -(n -1), (4)
which shows that the number of eigenvalues grows exponentially with respect to the dimension n. We have π 1 = 1, π 2 = 3 while π 3 is either 9 or 10. We note also for example that π 20 ≥ 1 572 861 and that obtaining an explicit formula for π n or improving the lower and the upper bounds in ( 4) are open questions. The lower bound in ( 4) is not sharp. In fact, the matrix A 2 of order 4 (given in Sect. 6) has 23 Paretoeigenvalues, which means that π 4 ≥ 23, while the lower bound in (4) gives π 4 ≥ 21.

In many applications, the Pareto-eigenvalue λ is assumed to be > 0 (see for example [START_REF] Judice | The eigenvalue complementarity problem[END_REF][START_REF] Júdice | On the solution of the symmetric eigenvalue complementarity problem by the spectral projected gradient algorithm[END_REF][START_REF] Júdice | On the asymmetric eigenvalue complementarity problem[END_REF]). In the rest of the paper, we will assume that λ > 0. By shifting the matrix A, it is possible to find nonpositive Pareto-eigenvalues.

The Semismooth Newton Method SNM

We start this section by reminding the reader some basic tools from nonsmooth analysis, and some basic properties of the Semismooth Newton Method (SNM). For more details, we refer to [START_REF] Qi | Convergence analysis of some algorithms for solving nonsmooth equations[END_REF][START_REF] Qi | A nonsmooth version of Newton's method[END_REF][START_REF] Ulbrich | In: Semismooth Newton Methods for Variational Inequalities and Constrained Optimization Problems in Function Spaces[END_REF] and references therein.

Consider a mapping Φ : R n → R n supposed to be locally Lipschitz. The B-subdifferential of Φ at a point z ∈ R n is defined by

∂ B Φ(z) = M ∈ M n (R) : ∃(z k ) ⊂ D Φ : z k → z and lim k→+∞ ∇Φ(z k ) = M ,
where D Φ is the set of differentiability points of Φ. The Clarke generalized Jacobian [START_REF] Clarke | Optimization and Nonsmooth Analysis[END_REF] of Φ is given by

∂Φ(z) = co ∂ B Φ(z), (5)
where "co" stands for the convex hull of the set ∂ B Φ(z).

The function Φ is said to be semismooth [START_REF] Mifflin | Semismooth and semiconvex functions in constrained optimization[END_REF][START_REF] Qi | A nonsmooth version of Newton's method[END_REF] at z ∈ R n if it is locally Lipschitz around z, directionally differentiable at z and satisfies the following condition sup

M∈∂Φ(z+h) Φ(z + h) -Φ(z) -Mz = o h . ( 6
)
If we replace o( h ) by O( h 2 ) in (6), then Φ is said to be strongly semismooth at z.

We note that semismooth functions are between Lipschitz continuous functions and continuously differentiable functions. For example, the Euclidean norm and piecewise differentiable functions are semismooth. Therefore, the functions

φ : R 2 → R, (x 1 , x 2 ) → φ max (x 1 , x 2 ) = max(x 1 , x 2 ), φ min (x 1 , x 2 ) = min(x 1 , x 2 ) and φ FB (x 1 , x 2 ) = x 1 + x 2 -x 2
1 + x 2 2 are semismooth. Given a locally Lipschitz function Φ : R n → R n , we consider the following nonlinear system of equations

Φ(x) = 0 R n . ( 7
)
We note that many problems in engineering and sciences can be formulated as [START_REF] De Luca | A semismooth equation approach to the solution of nonlinear complementarity problems[END_REF].

Next, we present the standard algorithm of the Semismooth Newton Method SNM, which will be used later.

Algorithm SNM

Initialization. Choose an initial point z 0 and set k = 0. Iteration. 1 One has a current point z k . Choose M k ∈ ∂Φ(z k ) and compute h k by solving the linear system

M k h k = -Φ z k . ( 8
)
Then, set

z k+1 = z k + h k and k = k + 1.
To ensure that the linear system (8) admits a unique solution, the matrix M k must be nonsingular. Furthermore, the following theorem gives some conditions to ensure the convergence of the above algorithm [START_REF] Qi | A nonsmooth version of Newton's method[END_REF].

Theorem 1 Let z be a zero of the function Φ. Suppose the following (i) Φ is semismooth (resp. strongly semismooth) at z; (ii) all matrices in ∂Φ(z) are nonsingular.

Then, there exists a neighborhood V of z such that the SNM initialized at any z 0 ∈ V generates a sequence (z k ) k∈N that converges superlinearly (resp. quadratically) to z.

Remark 2 Theoretical convergence results of the semismooth Newton method are ensured under the crucial assumption (ii) that all generalized Jacobian are nonsingular at the solution z. From a computational point of view, the rate of convergence of the SNM will also depend on the conditioning number of the matrices in ∂Φ(z) when z is around z.

As a first step towards a reformulation of the EiCP defined in (1), as a system of equations, we write

x ≥ 0, y≥ 0, x, y = 0, λx-Ax -y = 0, 1 n , x -1 = 0, (9)
1 As stopping criteria, we use Φ(z k ) < 10 -8 .

where 1 n is a vector of ones. The last condition in ( 9) is used to ensure that x is a nonzero vector. One can also use some other normalization conditions as x 2 2 = 1. Using the NCP-functions, it is possible to write the first three conditions in [START_REF] Dirkse | The path solver: a non-monotone stabilization scheme for mixed complementarity problems[END_REF] in the form U φ (x, y) = 0 R n where U φ : R n × R n → R n is the vector function defined by

U φ (x, y) = ⎡ ⎢ ⎣ φ(x 1 , y 1 ) . . . φ(x n , y n ) ⎤ ⎥ ⎦ ,
and φ : R 2 → R is a NCP-function, i.e. φ(a, b) = 0 ⇐⇒ a ≥ 0, b ≥ 0, ab = 0. ( 10
)
More precisely, we are led to solve the following system of 2n + 1 equations

U φ (x, y) = 0 R n , ( 11
) λx -Ax -y = 0 R n , ( 12
)
1 n , x -1 = 0. ( 13
)
In this paper, we focus on the NCP-functions, the most widely used in the literature, namely the Fischer-Burmeister (FB) function [START_REF] Fischer | A special Newton-type optimization method[END_REF] and the natural residual function [START_REF] Harker | Finite-dimensional variational inequality and nonlinear complementarity problem: a survey of theory, algorithms and applications[END_REF][START_REF] Pang | Newton's method for B-differentiable equations[END_REF], also called the min function, defined respectively by

φ FB (a, b) = a + b -a 2 + b 2 , ( 14
)
φ min (a, b) = min(a, b). ( 15
)
We note that other NCP-functions [START_REF] De Luca | A semismooth equation approach to the solution of nonlinear complementarity problems[END_REF][START_REF] Kanzow | New NCP-functions and their properties[END_REF][START_REF] Qi | Regular pseudo-smooth NCP and BVIP functions and globally and quadratically convergent generalized Newton methods for complementarity and variational inequality problems[END_REF] can be treated in the same way. Consider the vector function

Φ : R n × R n × R → R 2n+1 , z = (x, y, λ) -→ Φ(x, y, λ) defined by Φ(z) = Φ(x, y, λ) = ⎡ ⎣ U φ (x, y) λx -Ax -y 1 n , x -1 ⎤ ⎦ . ( 16
)
Therefore, the EiCP can be reformulated as a system of nonlinear and nonsmooth equations, i.e. (x, y, λ) is a solution of problem (1) if and only if

Φ(x, y, λ) = 0 R 2n+1 . ( 17
)

Lemma 2

The function Φ defined in ( 16) is semismooth with φ being one of the NCPfunctions φ min or φ FB . Moreover, its Clarke generalized Jacobian at z = (x, y, λ) is given by

∂Φ(z) = ⎧ ⎨ ⎩ ⎡ ⎣ E F 0 λI n -A -I n x 1 T n 0 0 ⎤ ⎦ : [E, F ] ∈ ∂U φ (x, y) ⎫ ⎬ ⎭ , ( 18
)
with I n representing the identity matrix of order n.

E = diag(a 1 , . . . , a n ), F = diag(b 1 , . . . , b n ) ∈ M n (R)
are diagonal matrices whose ith diagonal element is given by the following two cases:

case 1: if φ = φ FB , then, for i = 1, . . . , n, ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ a i = 1 - x i (x i , y i ) b i = 1 - y i (x i , y i ) if (x i , y i ) = (0, 0), ( 19
)
and

a i = 1 -ξ i b i = 1 -ρ i if (x i , y i ) = (0, 0), ( 20
)
for every

(ξ i , ρ i ) ∈ R 2 such that (ξ i , ρ i ) ≤ 1. case 2: if φ = φ min , then, for i = 1, . . . , n, ⎧ ⎪ ⎨ ⎪ ⎩ a i = 0, b i = 1 if x i > y i , a i = 1, b i = 0 if x i < y i , a i = α i , b i = 1 -α i , if x i = y i such that α i ∈ [0, 1]. ( 21
)
Proof The proof of ( 18) is given in [START_REF] Adly | A nonsmooth algorithm for cone-constrained eigenvalue problems[END_REF] while the proofs of ( 19), ( 20) and ( 21) can be found in [START_REF] Pang | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF].

The Lattice Projection Method (LPM)

In this section, we introduce a new formulation of the EiCP which leads us to solve a nonlinear and nonsmooth system of (2n + 1) equations involving the same number of variables. The originality of this formulation, in comparison with the existing literature, is that it is not based on the complementarity approach, i.e. we do not use any NCP-functions. We transform the EiCP to find the eigenvalues of a nonlinear problem. The following lemma is in this sense. Without loss of generality, we will assume in the sequel that λ > 0. If the matrix A has a nonpositive eigenvalue λ < 0, we set à = A + μI n with μ > 0 and large enough. It is easy to see, in this case that the spectrum σ ( Ã) ⊂ R * + .

Lemma 3

The EiCP defined in (1) with λ > 0 is equivalent to find the roots of the following nonlinear and nonsmooth function f :

R n × R * + → R n defined by (x, λ) -→ f (x, λ) = (Ax) + -λx, ( 22
)
where (Ax) + = max(Ax, 0) is the positive componentwise.

Proof Suppose λ > 0 and x ∈ R n \ {0}. EiCP defined in (1)

0 ≤ x ⊥ (λx -Ax) ≥ 0 is equivalent to the following problem min(x, λx -Ax) = 0,
where the min is taken componentwise. Since λ > 0, we have the following

min(x, λx -Ax) = 0 ⇐⇒ min(λx, λx -Ax) = 0 ⇐⇒ max(-λx, Ax -λx) = 0 ⇐⇒ max(Ax, 0) = λx,
which completes the proof.

Remark 3 Lemma 3 showed that 0 ≤ x ⊥ (λx -Ax) ≥ 0 ⇐⇒ (Ax) + = λx,
which means that the EiCP is equivalent to solve the nonlinear eigenvalue problem

(P R n + • A)(x) = λx,
where P R n + is the projection operator over R n + .

In the rest of this section, we are led to solve the following system of 2n

+ 1 equations ⎧ ⎪ ⎨ ⎪ ⎩ ỹ+ -λx = 0 R n , Ax -ỹ = 0 R n , 1 n , x -1 = 0.
We use the normalization condition 1 n , x = 1 to ensure that x is a nonzero eigenvector.

Let us consider the following function

Φ LPM : R n × R n × R → R 2n+1 defined by Φ LPM (x, ỹ, λ) = ⎡ ⎣ ỹ+ -λx Ax -ỹ 1 n , x -1 ⎤ ⎦ . ( 23
)
It is clear that solving EiCP is equivalent to solving the nonlinear system

Φ LPM (x, ỹ, λ) = 0 R 2n+1 . ( 24
)
We will use the SNM Algorithm presented in Sect. 3 to solve the generalized equation [START_REF] Martins | Dynamic stability of finite dimensional linearly elastic systems with unilateral contact and Coulomb friction[END_REF]. This method will be called Lattice Projection Method (LPM). In Sects. 6 and 7, the LPM will be compared to SNM FB and SNM min .

The following lemma gives a brief description of the representation of the Clarke generalized Jacobian of Φ LPM .

Lemma 4

The function Φ LPM , defined in [START_REF] Kettenring | Canonical analysis of several sets of variables[END_REF], is semismooth. Moreover, its Clarke generalized Jacobian at z = (x, ỹ, λ) is given by the following block representation

∂Φ LPM (z) = ⎧ ⎨ ⎩ ⎡ ⎣ -λI n F -x A -I n 0 1 T n 0 0 ⎤ ⎦ : F ∈ ∂ (•) + ( ỹ) ⎫ ⎬ ⎭ , ( 25
)
where F = diag(u 1 , . . . , u n ) ∈ R n×n is a diagonal matrix whose ith diagonal element is given by

u i = 1 + Sign( ỹi ) 2 = ⎧ ⎪ ⎨ ⎪ ⎩ 1 if ỹi > 0, [0, 1] if ỹi = 0, 0 if ỹi < 0. ( 26
)
Proof Follows immediately by using Theorem (Chain Rule II) in [START_REF] Clarke | Optimization and Nonsmooth Analysis[END_REF].

Nonsingularity conditions

In this section, we study under which conditions all the elements of the Clarke generalized Jacobians ∂Φ(z * ) defined in [START_REF] Judice | The eigenvalue complementarity problem[END_REF] and ∂Φ LPM (z * ) defined in [START_REF] Martins | Stability of finite-dimensional nonlinear elastic systems with unilateral contact and friction[END_REF] at solutions z * and z * of ( 17) and ( 24), respectively, are nonsingular. We begin first by fixing some notations that will be used later. The index set {1, 2, . . . , n} will be abbreviated by the capital letter

I . If M = (m ij ) ∈ M n (R) is a matrix and J , K ⊆ I , then M J K denotes the submatrix in M |J |,|K| (R) with elements m ij , i ∈ J , j ∈ K.
We recall that a matrix M ∈ M n (R) is a P-matrix if all its principal minors are strictly positive. Furthermore, M is a P-matrix if and only if for every nonzero vector x ∈ R n \ {0}, there exist an index i ∈ I such that x i = 0 and x i [Mx] i > 0.

Let us now introduce the function

F : R n × R n × R → R 2n+1 defined by z = (x, y, λ) -→ F (z) = ⎡ ⎣ y λx -Ax -y 1 n , x -1 ⎤ ⎦ . It is easy to see that F is of class C 1 .
Proposition 1 Any H ∈ ∂Φ(z), with φ being one of the NCP-functions φ min or φ FB , can be written in the form

H = D a + D b F (z), ( 27
)
where D a , D b ∈ M 2n+1 (R) are two positive semidefinite diagonal matrices defined by

D a = diag(a i ) 1≤i≤2n+1 and D b = diag(b i ) 1≤i≤2n+1 ,
whose ith diagonal element being one of the two cases mentioned above in [START_REF] Júdice | On the solution of the symmetric eigenvalue complementarity problem by the spectral projected gradient algorithm[END_REF], [START_REF] Júdice | On the asymmetric eigenvalue complementarity problem[END_REF] and [START_REF] Kanzow | A new class of semismooth Newton type methods for nonlinear complementarity problems[END_REF] for all i = 1, . . . , n and

a i = 0, b i = 1 for i = n + 1, . . . , 2n + 1.
Proof Let H ∈ ∂Φ(z). By Lemma 2, we have

H = ⎡ ⎣ E F 0 λI n -A -I n x 1 T n 0 0 ⎤ ⎦ .
On the other hand, we have

F (z) = ⎡ ⎣ 0 I n 0 λI n -A -I n x 1 T n 0 0 ⎤ ⎦ . ( 28
)
Hence, from these equalities and using [21, Proposition 2.6], the statement of the proposition follows easily and we get

H = D a + D b F (z).
To deal with the nonsingularity of all Clarke generalized Jacobian in ∂Φ(z * ) at a solution z * = (x * , y * , λ * ), let us introduce some notations.

The function

G : R n × R → R n+1 is defined by w = (x, λ) -→ G(w) = λx -Ax 1 n , x -1 . ( 29
)
We set

w * = x * , λ * ,
and

J = Δ a + Δ b G w * , ( 30
)
where

Δ a = diag(a 1 , . . . , a n , 0 R ), Δ b = diag(b 1 , . . . , b n , 1) two diagonal matrices of M n+1 (R) with (x * i , y * i ) = (x * i , (λ * x * -Ax * ) i ) at a solution z * = (x * , y * , λ * ) and G (w * ) is given by G w * = λI * n -A x * 1 T n 0 . ( 31
)
We note that

λ * Δ a = Δ λ * a ,
where

λ * a = λ * a 1 , λ * a 2 , . . . , λ * a n .
We set

J = Δ λ * a + Δ b G w * . ( 32
)
The following standard index sets will be used:

α := i ∈ {1, . . . , n, n + 1} : w * i > 0 = G i w * , β := i ∈ {1, . . . , n, n + 1} : w * i = 0 = G i w * , γ := i ∈ {1, . . . , n, n + 1} : w * i = 0 < G i w * .
Proposition 2 Let z * = (x * , y * , λ * ) be a fixed root of Φ defined in [START_REF] Hotelling | The most predictable criterion[END_REF] with w * = (x * , λ * ). If the submatrix G (w * ) αα is nonsingular and its Schur-complement

G w * ββ -G w * βα G w * -1 αα G w *
αβ is a P-matrix, then the Jacobians J and J defined respectively in [START_REF] Pang | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF] and (32) are nonsingular matrices.

Proof It is easy to see that the EiCP is equivalent to the following system

find w ∈ R n+1 \ {0}, such that w ≥ 0, G(w) ≥ 0, w, G(w) = 0, ( 33
)
where G is given by ( 29).

The desired result follows immediately from Lemma 5.3 in [START_REF] De Luca | A semismooth equation approach to the solution of nonlinear complementarity problems[END_REF] or Theorem 2.8 in [START_REF] Kanzow | A new class of semismooth Newton type methods for nonlinear complementarity problems[END_REF].

Theorem 2 Assume that z * = (x * , y * , λ * ) is a fixed root of Φ and set w * = (x * , λ * ).
Then the Jacobian J defined in [START_REF] Pang | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF] is nonsingular if and only if all elements H in ∂Φ(z * ) are nonsingular.

Proof Let H ∈ ∂Φ(z * ) and let X = (p, q, r) ∈ R n × R n × R be an arbitrary vector such that

H X = 0 R 2n+1 . ( 34
)
The previous equation ( 34) can be written explicitly as

⎡ ⎣ E F 0 λ * I n -A -I n x * 1 T n 0 0 ⎤ ⎦ ⎡ ⎣ p q r ⎤ ⎦ = ⎡ ⎣ 0 R n 0 R n 0 R ⎤ ⎦ ,
where E and F are defined in Lemma 2. Therefore,

Ep + F q = 0 R n , ( 35
)
λ * I n -A p -q + rx * = 0 R n , ( 36
)
1 n , p = 0. (37)
Using [START_REF] Queiroz | The symmetric eigenvalue complementarity problem[END_REF], we have

q = λ * I n -A p + rx * . ( 38
)
Consequently, [START_REF] Qi | Regular pseudo-smooth NCP and BVIP functions and globally and quadratically convergent generalized Newton methods for complementarity and variational inequality problems[END_REF] becomes

Ep + F λ * I n -A p + rF x * = 0 R n . ( 39
)
Therefore,

Ep + F λ * I n -A p + rF x * = 0 R n , 1 n , p = 0,
which can be rewritten in the following matricial form

E + F (λ * I n -A) F x * 1 T n 0 p r = 0 R n 0 R . ( 40
)
By setting Y = (p, r) ∈ R n × R, it is clear that ( 40) is equivalent to

J Y = 0 R n+1 . ( 41
)
Hence,

H X = 0 R 2n+1 ⇐⇒ J Y = 0 R n+1 , q = λ * I n -A p + rx * .
The desired conclusion of Theorem 2 follows immediately from the last equivalence.

The following theorem gives nonsingularity conditions of the generalized Jacobian of Φ LPM .

Theorem 3 Assume that z * = (x * , ỹ * , λ * ) is a fixed root of Φ LPM defined in [START_REF] Kettenring | Canonical analysis of several sets of variables[END_REF] and set w * = (x * , λ * ). Then the Jacobian J defined in [START_REF] Pinto Da Costa | Cone-constrained eigenvalue problems: theory and algorithms[END_REF] is nonsingular if and only if all elements in ∂Φ LPM (z * ) are nonsingular.

Proof Let H ∈ ∂Φ LPM (z * ) and let X = (p, q, r) ∈ R n × R n × R be an arbitrary vector such that

H X = 0 R 2n+1 . ( 42
)
The previous equation ( 42) can be written explicitly as

⎡ ⎣ -λ * I n F -x * A -I n 0 1 T n 0 0 ⎤ ⎦ ⎡ ⎣ p q r ⎤ ⎦ = ⎡ ⎣ 0 R n 0 R n 0 R ⎤ ⎦ ,
where F is defined in Lemma 4. Therefore,

-λ * p + F q -rx * = 0 R n , ( 43
)
Ap -q = 0 R n , ( 44
)
1 n , p = 0. ( 45
)
Using (44), we have

q = Ap. ( 46
)
Consequently, [START_REF] Ulbrich | In: Semismooth Newton Methods for Variational Inequalities and Constrained Optimization Problems in Function Spaces[END_REF] becomes

λ * I n -F A p + rx * = 0 R n . ( 47
)
Therefore,

λ * I n -F A p + rx * = 0 R n , 1 n , p = 0. ( 48
)
Notice that

E + F = I n where E = diag(a 1 , . . . , a n ), F = diag(b 1 , . . . , b n ) ∈ M n (R)
are diagonal matrices whose ith diagonal element is given by ( 21) and φ min being the min function defined in [START_REF] Horst | Relations among m sets of measures[END_REF]. Therefore, (48) becomes

λ * Ep + λ * F -F A p + rx * = 0 R n , 1 n , p = 0. ( 49
)
Since y * = λ * x * -Ax * , then we get ỹ * = Ax * = λ * x *y * . Hence, F = diag(u 1 , . . . , u n ) whose ith diagonal element is given by

u i = ⎧ ⎪ ⎨ ⎪ ⎩ 1 ifλ * x * i > y * i , [0, 1] if λ * x * i = y * i , 0 ifλ * x * i < y * i . (50)
Recall that F = diag(b 1 , . . . , b n ) is defined by

b i = ⎧ ⎪ ⎨ ⎪ ⎩ 1 ifx * i > y * i , [0, 1] if x * i = y * i , 0 ifx * i < y * i .
(51)

Using the complementarity approach 0 ≤ x * ⊥ y * ≥ 0, and the fact that λ * > 0, we get 0 ≤ λ * x * ⊥ y * ≥ 0 and x * i = 0 or y * i = 0.

Therefore,

F = F.
We note also that we have obviously

F x * = x * (since E + F = I n ).
Consequently, (49) becomes

λ * Ep + F λ * I n -A p + rx * = 0 R n , 1 n , p = 0, (52)
which can be rewritten in the following matricial form

Ẽ + F (λ * I n -A) F x * 1 T n 0 p r = 0 R n 0 R , ( 53
)
where

Ẽ = λ * E. By setting Y = (p, r) ∈ R n × R, it is clear that (53) is equivalent to J Y = 0 R n+1 . ( 54
)
Hence,

H X = 0 R 2n+1 ⇐⇒ J Y = 0 R n+1 , q = Ap.
The desired conclusion of Theorem 3 follows immediately from the last equivalence.

In order to illustrate the above results, we consider the following examples.

Example 1 Consider the following matrix

A = 8 -1 3

.

Using Lemma 1, it is easy to see that λ * = 4 is a solution of the EiCP with x * = (0, 1) T as an associated eigenvector. We set z * = (0, 1, 1, 0, 4) T , z * = (0, 1, -1, 4, 4), and we show that z * (respectively z * ) is a "good solution" of the nonlinear equation Φ(z) = 0 defined in (17) (respectively Φ LPM (z) = 0 defined in ( 24)).

A direct application of Proposition 2 yields

G w * = ⎡ ⎣ -4 1 0 -3 0 1 1 1 0 ⎤ ⎦ and G w * αα = 0 1 1 0 (which is a nonsingular matrix),
with w * = (0, 1, 4) T , α = {2, 3}, β = ∅ and γ = {1}. Hence, the first condition of Proposition 2 is satisfied and so is the second, since β is empty. Consequently, J and J are nonsingular matrices and so are H and H (by Theorems 2 and 3).

On the other hand, a direct calculation of the Jacobian matrix H ∈ ∂Φ(z * ) shows that it is nonsingular and is given by

H = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 1 0 0 0 0 0 0 0 1 0 -4 1 -1 0 0 -3 0 0 -1 1 1 1 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦
, with φ being one of the NCP-functions φ min or φ FB . Moreover, we have

F = F = 0 0 0 1 ,
where ỹ * = Ax * = (-1, 4) T . A direct calculation of the Jacobian matrix H ∈ ∂Φ LPM (z * ) shows that it is nonsingular and is given by

H = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ -4 0 0 0 0 0 -4 0 1 -1 8 -1 -1 0 0 3 4 0 -1 0 1 1 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ . Example 2 Let B = ⎡ ⎣ 15 -2 0 -11 9 -1 6 -10 23 ⎤ ⎦ .
Lemma 1 ensures that λ * = 23 is a solution of EiCP associated to B with x * = (0, 0, 1) T as an eigenvector. In this case, y * = λ * x * -Bx * = (0, 1, 0) T and ỹ * = Bx * = (0, -1, 23) T . We set z * = (0, 0, 1, 0, 1, 0, 23) T , z * = (0, 0, 1, 0, -1, 23, 23) T , and we show that z * (respectively z * ) is a "good solution" of the nonlinear equation

Φ(z) = 0 (respectively Φ LPM (z) = 0). A direct application of Proposition 2 yields G w * = ⎡ ⎢ ⎢ ⎣ 8
G w * ββ -G w * βα G w * -1 αα G w * αβ = 8 -0 0 0 1 1 0 -6 1 = 8 > 0,
is obviously a P-matrix. Hence, J and J are nonsingular and so are H and H by Theorems 2 and 3. On the other hand, a direct calculation of the Jacobian matrix H , with φ min being the min function defined in [START_REF] Horst | Relations among m sets of measures[END_REF], yields

H = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ δ 0 0 1-δ 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 8 2 0 -1 0 0 0 11 14 1 0 -1 0 0 -6 10 0 0 0 -1 1 1 1 1 0 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , where δ ∈ [0, 1]. The determinant of H is given by det(H ) = 8 -7δ.
We note that H is nonsingular since δ ∈ [0, 1].

Similarly, the determinant of the Jacobian matrix

H H = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ -23 0 0 δ 0 0 0 0 -23 0 0 0 0 0 0 0 -23 0 0 1 -1 15 -2 0 -1 0 0 0 -11 9 -1 0 -1 0 0 6 -10 23 0 0 -1 0 1 1 1 0 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , is given by det( H ) = -23(23 -15δ).
Hence, H is nonsingular since δ ∈ [0, 1].

Example 3 Let C = 1 -1 0 1 .
Analogously to the two examples above, using Lemma 1, λ * = 1 is a solution of EiCP associated to C with x * = (1, 0) T as an eigenvector, y * = λ * x * -Cx * = (0, 0) T and ỹ * = Cx * = (1, 0) T . We show that z * = (1, 0, 0, 0, 1) T and z * = (1, 0, 1, 0, 1) T are "bad solutions". A direct calculation of H ∈ ∂Φ(z * ) shows that it is not always nonsingular. If we take the NCP-function φ min defined in [START_REF] Horst | Relations among m sets of measures[END_REF], we get

H = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 1 0 0 0 δ 0 1-δ 0 0 1 -1 0 1 0 0 0 -1 0 1 1 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ , where δ ∈ [0, 1]. We have det(H ) = -δ,
which means that H is nonsingular if and only if δ ∈ ]0, 1].

In the same way, since

det( H ) = 1 -δ,
where δ ∈ [0, 1], and

H = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ -1 0 1 0 1 0 -1 0 δ 0 1 -1 -1 0 0 0 1 0 -1 0 1 1 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦
, then the Jacobian matrix H is nonsingular if and only if δ ∈ [0, 1[. On the other hand, using Proposition 2, we get

G w * = ⎡ ⎣ 0 1 1 0 0 0 1 1 0 ⎤ ⎦ and G w * αα = 0 1 1 0
(which is a nonsingular matrix), with α = {1, 3}, β = {2} and γ = ∅. However, since the Schur-complement of G αα (w * ) is equal to zero, it is not a P-matrix. Then, the second condition of Proposition 2 is not satisfied.

Testing on matrices of order 3, 4 and 5

The first numerical experiment will be given by taking matrices of order 3, 4 and 5 that are known to have 9, 23 and 57 Pareto-eigenvalues, respectively. We compare LPM, SNM FB and SNM min by computing the number of iterations, the CPU-time and the number of failures.

A 1 = ⎡ ⎣ 5 -8 2 -4 9 1 -6 -1 13 ⎤ ⎦ , A 2 = ⎡ ⎢ ⎢ ⎣ 132 -
Remark 4 In our case, a failure is declared if the Jacobian matrix is ill-conditioned or in the SNM, the number of iterations exceeded 100.

The comparison results are summarized in Table 1 where "Iter" denotes the average iterative number,"Time" denotes the average CPU-time and "Failure (%)" represents the percentage of the number of failures to find a solution for (EiCP).

We have used 10 5 random initial points, in order to obtain the 57 Paretoeigenvalues of A 3 simultaneously by the three solvers. We note that LPM finds all eigenvalues of A 3 by using only 10 3 random initial points.

Conclusion

As a first numerical experiment, it is clear that the LPM captures our attention by comparing the percentage of failures of each solver since it has the least percentage. We can also observe that the number of iterations of the SNM FB is high compared to the two others.

In the next section, we confirm our result by enlarging the set of tested examples. Moreover, we will compare the CPU time of the three solvers, since for the moment, LPM and SNM min have almost the same computing time.

Numerical results

In this section, we compare the three solvers that have been defined in Sects. 3 and4, namely LPM, SNM FB and SNM min . In order, to complete this experience, we choose the performance profiles developed by E.D. Dolan and J.J. Moré [6, 10] as a tool for comparing the solvers. The performance profiles give, for each t ∈ R, the proportion ρ s (t) of test problems on which each solver under comparison has a performance within the factor t of the best possible ratio.

Computing time, number of iterations, number of failures and maximum number of eigenvalues found by each solver are used as performance measures to compare these algorithms.

Due to the absence of library dedicated to EiCP, we have chosen a set P of 100 matrices, most of them are taken from the Matrix Market [START_REF] Boisvert | Matrix market: a web resource for test matrix collections[END_REF], with size less than 350. Let S be the set of the three solvers that will be compared and n s the number of solvers. The performance ratio is defined by

r p,s = t p,s min{t p,s : s ∈ S} ,
where p ∈ P , s ∈ S, and t p,s is either

• the average CPU time required to solve problem p by solver s corresponding to Fig. 1, or • the maximum number of solutions corresponding to Fig. 2, or • the number of failures (in the sense of Remark 4) corresponding to Fig. 3, or • the average iterative number corresponding to Fig. 4.

The performance of the solver s ∈ S is defined by where, n p is the number of problems, and t is a real factor. In this case, ρ s (t) is the probability, for solver s ∈ S, that the performance ratio r p,s is within a factor t ∈ R of the best possible ratio. For more details, we refer to [START_REF] Dolan | Benchmarking optimization software with performance profiles[END_REF]. Thanks to the performance profile plot, one can compare how well a solver can perform relative to the other by taking in consideration the four criteria listed above. The value of ρ s [START_REF] Adly | A nonsmooth algorithm for cone-constrained eigenvalue problems[END_REF] gives the probability that the solver s will win over the two others, while ρ s (t), for large value of t , measures the robustness of the solver s.

ρ s (t) = 1 n p size {p ∈ P : r p,s ≤ t},
The numerical experiments are carried out in a Powerbook Mac OS 10.6.8 with a processor 2.33 GHz Intel Core 2 Duo and 2 Go memory. All the program codes are written and executed in Matlab 7.7.

Figures 1,2,3, and 4 show the performance profiles of each solver and a detailed description of the tested problems is given.

Figure 1 presents the performance profiles of the three solvers corresponding to the average computing time. We observe that the LPM is efficient and robust. In fact, in the interval [0, 1], LPM can solve 99 % of the problems, while SNM FB and SNM min do not reach 40 % and require a long running time. SNM FB has a lower number of wins than either LPM and SNM min . However, for t ≥ 2, the three solvers confirm their robustness.

Figure 1 also indicates that relative to the computing time, with the same initial points and under the same stopping criterion, LPM is the fastest solver, followed respectively by SNM min and SNM FB .

In Fig. 2, the performance ratio represents the maximum number of solutions found by each solver. We remark that the LPM is the winner, it is able to solve 99 % of the problems and detects the maximum number of Pareto-eigenvalues. SNM FB has a lower number of wins than either LPM or SNM min but its performance increases with the growth of the factor t .

The LPM captures our attention in Fig. 3, with its ability to solve over 99 % of problems. Performance ratio represents here the number of failures to solve each problem by each solver (see Remark 4). Another point of interest is that the number of failures giving by SNM min is important than the other solvers. We conclude again that the LPM is the winner followed respectively by SNM FB and SNM min .

In Fig. 4, performance ratio displays the average iterative number for each solver. It is clear that SNM min has the most wins, followed by LPM and SNM FB . We observe that the performance of LPM becomes interesting when t ≥ 2.

We conclude that LPM confirms again its efficiency, robustness and its ability to solve EiCP. By comparing the three solvers on a set of problems with respect to the computing time, number of iterations, number of failures and maximum number of found eigenvalues, LPM performs better than the two other solvers.

Applications to Pareto bivariate Eigenvalue Complementarity Problems (bi-EiCP)

We define, in this section, the so-called, bivariate Eigenvalue Complementarity Problems (bi-EiCP for short), in the case of Pareto cones. Surprisingly, the numerical simulations for this subject have never been studied. We propose to generalize the LPM, defined in Sect. 4, and to use the same tools mentioned in the previous section as performance measures to compare our method with SNM FB and SNM min for solving bi-EiCP. (

) 55
The problem (55) can be rewritten in the equivalent form

Ez = Λz, z i = 1, i = 1, 2, (56)
where

E = A B C D , Λ:= diag{λI n , μI m } and z = z 1 z 2 = x y .
Bi-EiP arises in many fields of applied mathematics and engineering [START_REF] Kettenring | Canonical analysis of several sets of variables[END_REF][START_REF] Tanaka | Some generalized methods of optimal scaling and their asymptotic theories: the case of multiple responses-multiple factors[END_REF]. They are necessary for the understanding of the constrained bi-EiCP. Bi-EiP was first studied by Hotelling in a brief paper [START_REF] Hotelling | The most predictable criterion[END_REF] in 1935. His definitive study of the problem appeared one year later in [START_REF] Hotelling | Relations between two sets of variates[END_REF] and still stands as a key reference in multivariate statistical literature. In 1961, an iterative method, proposed by Horst [START_REF] Horst | Relations among m sets of measures[END_REF], has been used for solving bi-EiP. Three decades later, Chu and Watterson [START_REF] Chu | On a multivariate eigenvalue problem: I. Algebric theory and a power method[END_REF] proved the convergence of Horst's iterative method. Among recent contributions of bi-EiP, we can cite [START_REF] Hanafi | Global optimality of the successive Maxbet algorithm[END_REF] and references therein. In this section, we propose the use of SNM to solve bi-EiCP.

Pareto bi-Eigenvalue Complementarity Problems For real matrices

A ∈ M n (R), B ∈ M n,m (R), C ∈ M m,n (R) and D ∈ M m (R), the bi-EiCP associated to the Pareto cones (R n + , R m + ) is written as follows: ⎧ ⎪ ⎨ ⎪ ⎩ find (λ, μ) ∈ R × R and (x, y) ∈ R n \ {0} × R m \ {0}, such that x ≥ 0, λx -Ax -By ≥ 0, x, λx -Ax -By = 0, y ≥ 0, μy -Cx -Dy ≥ 0, y, μy -Cx -Dy = 0. (57)
The couple (λ, μ) ∈ R × R is called the Pareto-bi-eigenvalue of the block structured matrix E.

In [START_REF] Amri | Spectral analysis of coupled linear complementarity problems[END_REF], A. Amri and A. Seeger generalize Lemma 1 to solve bi-EiCP. They characterized the Pareto-bi-eigenvalue (57). The following Lemma is in this sense.

Lemma 5 [START_REF] Amri | Spectral analysis of coupled linear complementarity problems[END_REF] Given a pair (λ, μ) ∈ R × R, the following statements are equivalent:

(i) Problem (57) has a nonzero solution (x, y) ∈ R n × R m .
(ii) There are index sets I ⊆ N = {1, . . . , n} and J ⊆ M = {1, . . . , m} such that

A II ξ + B IJ η = λξ C JI ξ + D JJ η = μη (58) has a solution (ξ, η) ∈ R |I |+|J | satisfying the interiority conditions ξ ∈ int R |I | + , η∈ int R |J | + (59)
and the binding conditions

k∈I a ik ξ k + l∈J b il η l ≤ 0, ∀ i ∈ N \ I, ( 60
) k∈I c jk ξ k + l∈J d jl η l ≤ 0, ∀ j ∈ M \ J. (61)
Furthermore, (x, y) can be constructed from (ξ, η) by setting

x i = ξ i if i ∈ I, 0 if i ∈ N \ I, and y j = η j if j ∈ J, 0 if j ∈ M \ J. ( 62
)
Without loss of generality, we will assume that λ, μ > 0. In the same manner as EiCP, as a first step towards a reformulation of bi-EiCP as a system of equations, we write

x ≥ 0, y ≥ 0, u ≥ 0, v ≥ 0, u T x = 0, v T y = 0, λx -Ax -By -u = 0, μy-Cx -Dy -v = 0, 1 n , x -1 = 0, 1 m , y -1 = 0.
Using NCP-functions, the first three conditions above can be rewritten in the succinct form U ϕ (x, u) = 0 and V ϕ (y, v) = 0, where U ϕ : R 2n → R n and V ϕ : R 2m → R m are respectively given by

U ϕ (x, u) = ⎡ ⎢ ⎣ ϕ(x 1 , u 1 ) . . . ϕ(x n , u n ) ⎤ ⎥ ⎦ , V ϕ (y, v) = ⎡ ⎢ ⎣ ϕ(y 1 , v 1 ) . . . ϕ(y m , v m ) ⎤ ⎥ ⎦ .
Here, ϕ is one of the NCP-functions ( 14) or [START_REF] Horst | Relations among m sets of measures[END_REF].

In short, we are led to solve the following system of d = 2(n + m + 1) equations

Φ(z) = 0 R d ,
where Φ : R d → R d is a locally Lipschitz function defined by

Φ(z) = Φ(x, u, y, v, λ, μ) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ U ϕ (x, u) V ϕ (y, v) λx -Ax -By -u μy -Cx -Dy -v 1 n , x -1 1 m , y -1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . ( 63
)
Lemma 6 The Clarke generalized Jacobian of Φ at z = (x, u, y, v, λ, μ) is given by (y, v), where E 1 , F 1 , E 2 and F 2 can be computed in the same way as in Lemma 2.

∂Φ(z) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ E 1 F 1 0 0 0 0 0 0 E 2 F 2 0 0 λI n -A -I n -B 0 x 0 -C 0 μI m -D -I m 0 y 1 T n 0 0 0 0 0 0 0 1 T m 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (64) with [E 1 , F 1 ] ∈ ∂U ϕ (x, u) and [E 2 , F 2 ] ∈ ∂V ϕ

Lattice Projection Method applied to bi-EiCP

In this section, we reformulate the bi-EiCP into a system of nonlinear equations. The following lemma is in this sense.

Lemma 7

Solving bi-EiCP is equivalent to find the roots of the following nonlinear and semismooth function Φ LPM : R d → R d defined by

Φ LPM (z) = Φ LPM (x, u, y, v, λ, μ) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ u + -λx Ax + By -u v + -μy Cx + Dy -v 1 n , x -1 1 n , y -1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . ( 65
)
The Clarke generalized Jacobian of Φ LPM at z = (x, u, y, v, λ, μ) is given by

∂Φ LPM (z) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ -λI n E 3 0 0 -x 0 A -I n B 0 0 0 0 0 -μI m F 3 0 -y C 0 D -I m 0 0 1 T n 0 0 0 0 0 0 1 T m 0 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , ( 66
)
where E 3 ∈ ∂(•) + (u) and F 3 ∈ ∂(•) + (v) (see Lemma 4).

Numerical experiments for solving bi-EiCP

Similarly to Sect. 7, we study here the following three solvers:

(i) SNM FB , where Φ : R d → R d is defined in (63) and ϕ is the Fischer-Burmeister function ( 14), (ii) SNM min , where Φ is defined in (63) and ϕ is the min-function ( 15), (iii) The generalization of the LPM, where Φ LPM is defined in (65). For the comparison, we may use computing time, number of failures, maximum number of solutions and number of iterations as performance measures. Due to the long time required for the simultaneous execution of solvers, we restrict our analysis to a set of 80 matrices for solving bi-EiCP.

In Fig. 5, LPM captures our attention since it has the most wins and it does not require a long running time to solve bi-EiCP. Relative to the computing time, LPM is the fastest solver, followed respectively by SNM FB and SNM min .

Figure 6 represents the performance profiles of the three solvers with respect to the number of failures. We observe that LPM is able to solve 99 % of problems, whereas SNM min has a lower number of wins.

In Fig. 7, LPM confirms its ability to detect the maximum number of solutions. We also note that the performance of SNM min is lower.

In Fig. 8, we compare the number of iterations corresponding to the common set of eigenvalues of each problem for the three solvers. With respect to this criterion, the performance of LPM and SNM min are very close, while SNM FB has a lower number of wins. We note that when t ≥ 0.5, the three solvers are comparable. Due to the difference of the performance between (LPM, SNM FB ) and SNM min , with respect to the maximum number of solutions (see Fig. 7), the performance given in Fig. 8 for the three solvers does not reach 86 %.

In Fig. 9, we represent only the performance profiles of LPM and SNM FB corresponding to the average iterative number. We note that both solvers solve 95 % of problems.

It is clear that, in this figure, the number of iterations of SNM FB is more important than LPM. SNM FB needs more iterations to solve each problem. LPM has the most wins, the robustness of the two solvers is confirmed when t ≥ 0.6.

Concluding remarks

In this paper, we proposed a new algorithm, namely the Lattice Projection Method (LPM), for solving Eigenvalue Complementarity Problems and bi-Eigenvalue Complementarity Problems. This method is compared to SNM min and SNM FB , two algorithms widely discussed in the literature for solving nonlinear complementarity problems. The performance profiles assessment, evaluated on more than 100 benchmark datasets, indicate the efficiency and the robustness of the LPM. It would be interesting to consider the globalization of the LPM and compare it to other solvers like PATH-solver [START_REF] Dirkse | A pathsearch damped Newton method for computing general equilibria[END_REF][START_REF] Dirkse | The path solver: a non-monotone stabilization scheme for mixed complementarity problems[END_REF]. This is out of the scope of this paper and will probably be the subject of a future project of research. Another forthcoming work will involve the extension of the proposed algorithm to the second order cone (also called the Lorentz cone).

  a nonsingular matrix), with w * = (0, 0, 1, 23) T , α = {3, 4}, β = {1} and γ = {2}. Moreover, the Schurcomplement of G (w * ) αα ,

Fig. 1 Fig. 2

 12 Fig. 1 Performance profiles where t p,s is the average computing time (Color figure online)

Fig. 3 Fig. 4

 34 Fig. 3 Performance profiles where t p,s is the number of failures (Color figure online)

8. 1

 1 The classical bi-Eigenvalue Problems (bi-EiP) Let A, B, C and D be real matrices of order n × n, n × m, m × n and m × m, respectively. The classical unconstrained bivariate Eigenvalue Problem is defined by(bi-EiP) ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ find a pair (λ, μ) ∈ R × R, (x, y) ∈ R n × R m , such thatAx + By = λx Cx + Dy = μy x = y = 1.

Fig. 5 Fig. 6

 56 Fig.[START_REF] Clarke | Optimization and Nonsmooth Analysis[END_REF] Performance profiles where t p,s is the average computing time (Color figure online)

Fig. 7 Fig. 8 Fig. 9

 789 Fig. 7 Performance profiles where t p,s is the maximum number of solutions (Color figure online)

Table 1

 1 Comparison of the three solvers on the matrices A 1 , A 2 and A 3

	106 18 74 24 101 ⎤ 81 -2 -44 195 -92 7 ⎥ ⎥ ⎦
	-21 -38	0	230

Acknowledgements

The first author wishes to thank Professor Jonathan Borwein for some stimulating discussions about the new method LPM and particularly Lemma 3.

The authors wish to thank also Professor Alberto Seeger for his careful reading of the first version of the manuscript.