Tecla De Luca
email: tdeluca@fub.it

Francisco Facchinei

Christian Kanzow
email: kanzow@math.uni-hamburg.de

A Theoretical and Numerical Comparison of Some Semismooth Algorithms for Complementarity Problems

Keywords: nonlinear complementarity problem, semismoothness, Newton's method, projected gradient method, large-scale problem

In this paper we introduce a general line search scheme which easily allows us to define and analyze known and new semismooth algorithms for the solution of nonlinear complementarity problems. We enucleate the basic assumptions that a search direction to be used in the general scheme has to enjoy in order to guarantee global convergence, local superlinear/quadratic convergence or finite convergence. We examine in detail several different semismooth algorithms and compare their theoretical features and their practical behavior on a set of large-scale problems.

Introduction

The main aim of this paper is twofold: On the one hand, we present a general scheme which allows us to easily analyze and compare, in a unified framework, the characteristics of several semismooth algorithms for the solution of the nonlinear complementarity problem x ≥ 0, F(x) ≥ 0, x T F(x) = 0, (NCP(F))

where F : IR n → IR n is a continuously differentiable function and all inequalities are taken componentwise. On the other hand, we give a numerical comparison of several realizations of our general framework. Although some of these realizations are known from the literature, some others lead to new algorithms. Interestingly, it turns out that these new semismooth algorithms have stronger theoretical properties than the known ones and that their numerical behaviour is at least comparable.

The nonlinear complementarity problem has important applications in operations research, economic equilibrium models and in the engineering sciences, see, e.g., [START_REF] Ferris | Engineering and economic applications of complementarity problems[END_REF][START_REF] Harker | Finite-dimensional variational inequality and nonlinear complementarity problems: A survey of theory, algorithms and applications[END_REF]. For this reason, there is a growing interest in finding efficient and robust algorithms for solving NCP(F). This reflects in an increasing number of proposals of solution schemes for NCP(F) in recent years. In these recent developments an important role has been played by the semismooth methods, i.e., by those methods that attempt to solve the complementarity problem by first reformulating it as a semismooth system of equations and then by applying a generalized Newton method to solve this system. We refer the reader to Subsection 2.1 for a short overview on semismooth Newton methods and to Subsection 2.3 for the presentation of some basic properties of two different reformulations which will be considered in this paper.

The first semismooth methods were based on the Fischer function [START_REF] Fischer | A special Newton-type optimization method[END_REF]:

ϕ F (a, b) := a 2 + b 2 -a -b.
The function ϕ F has the property that

ϕ F (a, b) = 0 ⇐⇒ a ≥ 0, b ≥ 0, ab = 0, (1)
so that the NCP(F) can be reformulated as a system of nonsmooth equations

F (x) = 0,
where F : IR n → IR n is defined by

F (x) := ⎛ ⎜ ⎜ ⎝ ϕ F (x 1 , F 1 (x)) . . . ϕ F (x n , F n (x)) ⎞ ⎟ ⎟ ⎠ .
The Fischer function, or some suitable modifications of it [START_REF] Chen | A penalized Fischer-Burmeister function: Theoretical investigation numerical results[END_REF][START_REF] Kanzow | A new class of semismooth Newton-type methods for nonlinear complementarity problems[END_REF][START_REF] Kanzow | New NCP-functions and their properties[END_REF][START_REF] Tseng | A new class of merit functions for the nonlinear complementarity problem[END_REF][START_REF] Sun | On NCP-functions[END_REF], has maintained a central role in the development of semismooth methods up to now; in fact it appears to possess many favourable properties both from the theoretical and numerical point of view.

In this paper we introduce a line search scheme based on the semismooth Newton method and on the Fischer function. This scheme encompasses several known methods and easily allows the development of new, interesting ones. More precisely, our general algorithmic scheme is a simple but useful generalization of the algorithms considered, e.g., in [START_REF] De Luca | A semismooth equation approach to the solution of nonlinear complementarity problems[END_REF][START_REF] Facchinei | A nonsmooth inexact Newton method for the solution of large-scale nonlinear complementarity problems[END_REF][START_REF] Facchinei | A new merit function for nonlinear complementarity problems and a related algorithm[END_REF]. It allows the use of a variety of different search directions; then it monitors convergence of the iterates by using the merit function

F (x) := 1 2 F (x) T F (x), (2)
i.e., the natural merit function of the operator F . We shall consider in detail several different specializations of the general scheme and compare the resulting algorithms both from the theoretical and numerical point of view. Another motivation for this work comes from the following observation: Almost all existing semismooth methods reformulate the complementarity problem as a system of equations (x) = 0 in order to obtain a locally fast convergent algorithm, and then use the corresponding merit function

(x) := 1 2 (x) T (x)
in order to globalize this method. This is a very natural approach which, usually, leads to very good algorithms. However, and this will also be discussed in more detail in Sections 2 and 3, sometimes one reformulation of the complementarity problem gives rise to a nicer merit function (e.g., weaker conditions for a stationary point to be a global minimum), whereas another reformulation gives the better search direction and therefore, in particular, the better local algorithm (e.g., finite termination or quadratic convergence under weaker assumptions). So it seems interesting to combine "the best" merit function with several schemes to compute a search direction. In fact, the probably most interesting algorithms among those studied in this paper are based on this combination idea. The organization of the paper is as follows: In Section 2, we restate some known properties of nonsmooth Newton methods, prove some preliminary results on certain regularity concepts and apply these results to two specific reformulations of the complementarity problem NCP(F). Our general algorithmic scheme is presented in Section 3. There, based on fairly general assumptions on the search directions computed within the algorithmic scheme, we also prove global and local convergence results. In Section 4, we consider six different realizations of our class of algorithms. We show that these realizations satisfy, under suitable conditions, all or most of the assumptions required for the search directions in the previous section. In Section 5, we then present a summary of our extensive numerical testing when using four of the six realizations of our general scheme. We conclude this paper with some final remarks in Section 6.

Some words about our notation: We say that F : IR n → IR m is a C 1 -function if F is continuously differentiable, and an LC 1 -function if F is differentiable with a locally Lipschitzian Jacobian. The Jacobian of a C 1 -function F : IR n → IR m at a point x ∈ IR n is denoted by F (x). If M ∈ IR n×n is any given matrix with elements m i j and J, K ⊆ {1, . . . , n} are two nonempty subsets, we write M JK for the submatrix consisting of the elements m i j with i ∈ J and j ∈ K . A similar notation will be used for subvectors.

Preliminaries

In this section we collect some material that will be used subsequently. In particular, in Subsection 2.1 we recall some known schemes for the solution of systems of semismooth equations along with their main convergence properties; these schemes will constitute the core step in the algorithms considered in this paper. In Subsection 2.2, in order to be able to better appreciate the theoretical differences between the algorithms we shall study, we elucidate the relations between b-regularity and R-regularity, the two assumptions most commonly used in analyzing convergence rates of algorithms for NCP(F). In particular, we point out the rather unexpected fact that, for P 0 -problems, R-regularity and b-regularity coincide. Finally, in Subsection 2.3 we apply these regularity concepts to prove nonsingularity results for two specific reformulations of the complementarity problem.

Local algorithms for semismooth systems

Let G : IR n → IR n be a locally Lipschitzian function. Then, by Rademacher's theorem, G is almost everywhere differentiable. If we denote by D G the set of points x ∈ IR n at which G is differentiable, we can define the B-subdifferential of G at x as

∂ B G(x) := H ∈ IR n×n | H = lim x k ∈D G ,x k →x G (x k ) ,
see [START_REF] Qi | Convergence analysis of some algorithms for solving nonsmooth equations[END_REF]. It is easy to see that this set is nonempty and compact. The convex hull of this set,

∂G(x) := conv∂ B G(x),
is the generalized Jacobian of G at x, see [START_REF] Clarke | Optimization and Nonsmooth Analysis[END_REF]. G is said to be semismooth at x [START_REF] Mifflin | Semismooth and semiconvex functions in constrained optimization[END_REF][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] if it is directionally differentiable at x and

H d -G (x; d) = o( d )
for any d → 0 and any H ∈ ∂G(x + d), where G (x; d) denotes the directional derivative of G at x in the direction d. Analogously, G is called strongly semismooth at x if

H d -G (x; d) = O( d 2 )
for any d → 0 and any

H ∈ ∂G(x + d). We call G (strongly) semismooth if it is (strongly) semismooth at any point x ∈ IR n . If G is a C 1 -function, then G is semismooth and ∂ B G(x) = ∂G(x) = {G (x)}
for every x ∈ IR n ; moreover, if G is an LC 1 -function, then G is strongly semismooth. This motivates the following generalization of the classical Newton method for the solution of the semismooth system of equations

G(x) = 0. ( 3
)
Generalized Newton Method

0. Choose x 0 ∈ IR n and set k := 0. 1. Select an element H k ∈ ∂ B G(x k ). Compute a solution d k ∈ IR n of the generalized Newton equation H k d = -G(x k ). (4) 2. Let x k+1 = x k + d k , k = k + 1 and go to Step 1.
Note that we do not incorporate any termination criterion in this and in the other algorithms presented in this section, since, in order to simplify the statement of their convergence properties, we always assume that they generate an infinite sequence of points. We call a solution x * ∈ IR n of (3) BD-regular if all elements H ∈ ∂ B G(x * ) are nonsingular, see [START_REF] Qi | Convergence analysis of some algorithms for solving nonsmooth equations[END_REF]. Now we can restate from [START_REF] Qi | Convergence analysis of some algorithms for solving nonsmooth equations[END_REF] the following local convergence result for the above algorithm.

Theorem 2.1. Assume that G : IR n → IR n is semismooth (strongly semismooth) and x * ∈ IR n is a BD-regular solution of (3). Then, for any x 0 ∈ IR n sufficiently close to x * , the above algorithm is well-defined and generates a sequence {x k } which converges to x * Q-superlinearly (Q-quadratically).

It is possible to consider variants of the Generalized Newton Method that are more suited to large-scale problems. The most obvious variant is an inexact version in which the Eq. ( 4) is only solved inexactly.

Generalized Inexact Newton Method

0. Choose x 0 ∈ IR n and set k := 0. 1. Select an element H k ∈ ∂ B G(x k ). Compute a solution d k ∈ IR n of the equation H k d = -G(x k ) + r k (5) for a suitable vector r k ∈ IR n . 2. Let x k+1 = x k + d k , k = k + 1 and go to Step 1.
The vector r k in ( 5) is called the residual and measures the inaccuracy with which the Newton Eq. ( 4) is solved. Note that usually, in actual computations, the vector r k is not fixed beforehands. Instead, an iterative method is used to solve the linear system (4), and this method is stopped when the norm of the residual r k is smaller than a prefixed accuracy. The following result can be established, see [START_REF] Facchinei | Inexact Newton methods for semismooth equations with applications to variational inequality problems[END_REF][START_REF] Martínez | Inexact Newton methods for solving nonsmooth equations[END_REF]. Theorem 2.2. Let G be semismooth and let x * be a BD-regular solution of the system G(x) = 0. Suppose that {x k } is a sequence generated by the Generalized Inexact Newton Method. Then the following statements hold: (a) There are numbers η > 0 and ε > 0 such that, if x 0x * ≤ ε and r k ≤ η G(x k ) for all k, the sequence {x k } is well-defined and converges Q-linearly to the solution x

H T k H k + σ k I d = -H T k G(x k ) + r k (6)
for a suitable vector r k ∈ IR n . 2. Let x k+1 = x k + d k , k = k + 1 and go to Step 1.

The vector r k is, again, the residual. The local properties of the Generalized Inexact Levenberg-Marquardt Method are also very similar to its smooth counterpart, and have been established in [START_REF] Facchinei | A nonsmooth inexact Newton method for the solution of large-scale nonlinear complementarity problems[END_REF].

Theorem 2.3. Let G : IR n → IR n be semismooth and let x * ∈ IR n be a BD-regular solution of G(x) = 0. Suppose that {x k } is a sequence generated by the Generalized Inexact Levenberg-Marquardt Method. Then the following statements hold: (a) There exist constants ε > 0, η > 0 and σ > 0 such that, if

x 0 -x * ≤ ε, r k ≤ η H T k G(x k ) and σ k ≤ σ for all k, then the sequence {x k } converges Q-linearly to x * . (b) If, in addition, r k = o( H T k G(x k ) ) and σ k → 0, then {x k } → x * Q-superlinearly. (c) If, in addition, G is strongly semismooth, and if r k = O( H T k G(x k ) 2 ) and σ k = O( H T k G(x k ) ), then {x k } → x * Q-quadratically.

On b-and R-regularity

Given a solution x * ∈ IR n of NCP(F), let us introduce the three index sets

α := {i ∈ I | x * i > 0 = F i (x * )}, β := {i ∈ I | x * i = 0 = F i (x * )}, γ := {i ∈ I | x * i = 0 < F i (x * )}.
The following two regularity concepts play an important role in the (theoretical and numerical) analysis of complementarity problems.

F (x * ) αα F (x * ) αβ F (x * ) βα F (x * ) ββ
is a P-matrix.

We recall that the above mentioned Schur-complement is defined by

F (x * ) ββ -F (x * ) βα F (x * ) -1 αα F (x * ) αβ ∈ IR |β|×|β| . ( 7
)
Obviously, there is no difference between band R-regularity at a nondegenerate solution of NCP(F). In general, however, it is known that b-regularity is a weaker assumption than R-regularity, see also Example 2.1 at the end of the next subsection. Here, we want to give a complete relationship between band R-regularity. To this end, we first prove a simple result.

Lemma 2.4. Let M ∈ IR n×n and σ ⊆ {1, . . . , n} be such that the submatrix M σ σ is nonsingular. Then the principal submatrices of the Schur-complement M/M σ σ are given by M σ ∪ν,σ ∪ν /M σ σ , where ν is any subset such that ∅ ⊆ ν ⊆ σ and σ := {1, . . . , n}\σ denotes the complementary subset of σ.

Proof: For any subset ν such that ∅ ⊆ ν ⊆ σ , we obtain

M M σ σ νν = M σ σ -M σ σ M -1 σ σ M σ σ νν = M σ ∩ν, σ ∩ν -M σ ∩ν,σ M -1 σ σ M σ, σ ∩ν = M νν -M νσ M -1 σ σ M σ ν = M σ ∪ν,σ ∪ν M σ σ ,
where the first and last equality follow directly from the definition of the Schur-complement and the third equation is a consequence of the fact that σ ∩ ν = ν. 2

We are now in the position to prove the main result of this subsection. (1) First assume that x * is an R-regular solution. Then the Schur-complement ( 7) is a P-matrix by definition. The fact that M α∪δ,α∪δ is nonsingular for all δ with ∅ ⊆ δ ⊆ β follows from the assumed R-regularity and [31, Lemma 1]. (2) Obviously, part (b) implies part (c) since any P-matrix is a P 0 -matrix.

(3) Now assume that x * is a b-regular solution of NCP(F) and that the Schur-complement S := M α∪β,α∪β /M αα is a P 0 -matrix. We have to show that this Schur-complement is actually a P-matrix, i.e., that det(S δδ ) > 0 for all subsets ∅ ⊆ δ ⊆ β. Since, by the assumed P 0 -property, we already know that det(S δδ ) ≥ 0, it suffices to show that det(S δδ ) = 0 for all ∅ ⊆ δ ⊆ β. In view of Lemma 2.4, we have

S δδ = M α∪δ,α∪δ M αα for any ∅ ⊆ δ ⊆ β.
Hence, from a well-known result on the determinant of a Schur-complement (see, e.g., [5, Proposition 2.3.5]), we get

det(S δδ ) = det M α∪δ,α∪δ M αα = det(M α∪δ,α∪δ ) det(M αα ) = 0,
where the last part follows from the fact that M α∪δ,α∪δ is nonsingular for any ∅ ⊆ δ ⊆ β by the assumed b-regularity of the solution x * . Hence x * is an R-regular solution of NCP(F). 2

The above result immediately implies the following corollary.

Corollary 2.6. Let x * ∈ IR n be a solution of NCP(F) such that the Jacobian F (x * ) is a P 0 -matrix. Then x * is an R-regular solution if and only if it is a b-regular solution.

Proof: The statement follows directly from Theorem 2.5 and the fact that the Schurcomplement of a nonsingular submatrix of a P 0 -matrix is again a

P 0 -matrix, see [3, Lemma 2.3]. 2
If F : IR n → IR n is a P 0 -function, then all Jacobian matrices are P 0 -matrices by [START_REF] Moré | On Pand S-functions and related classes of n-dimensional nonlinear mappings[END_REF]Theorem 5.8]. Hence there is no difference between Rand b-regular solutions for P 0function complementarity problems. In particular, there is no difference between these two regularity concepts for monotone problems. We finally note that the above results can also be obtained from more general ones about coherent orientation of the normal map [35].

The operators F and P

In this section we recall some basic (semi)smoothness properties of two operators which are often used in order to reformulate the nonlinear complementarity problem as a nonsmooth system of equations, namely the Minimum and the Fischer operator. The Fischer operator has already been defined in the introduction. The Minimum operator is similarly defined by

P (x) := ⎛ ⎜ ⎜ ⎝ min{x 1 , F 1 (x)} . . . min{x n , F n (x)} ⎞ ⎟ ⎟ ⎠ ;
here we use the subscript "P" for the Minimum operator in order to stress the importance of Pang's [START_REF] Pang | Newton's method for B-differentiable equations[END_REF][START_REF] Pang | A B-differentiable equation-based, globally and locally quadratically convergent algorithm for nonlinear programs, complementarity and variational inequality problems[END_REF] seminal work in this area. We believe that F (x) = 0 and P (x) = 0 are the most used reformulations of the complementarity problem as system of equations, and we shall use them extensively in this paper. In view of the results of the previous sections we are obviously interested in the semismoothness properties of these reformulations. The first result that can easily be established is that, if F is (strongly) semismooth, then also F and P are strongly semismooth, see [START_REF] Fischer | Solution of monotone complementarity problems with locally Lipschitzian functions[END_REF][START_REF] Kanzow | Semismooth Newton-type methods for the solution of nonlinear complementarity problems[END_REF] for details. In particular, we have the following Another point which deserves attention is the structure of the B-subdifferentials of F and P and the conditions which guarantee BD-regularity at a solution. The following result was shown in [START_REF] Facchinei | A new merit function for nonlinear complementarity problems and a related algorithm[END_REF].

Proposition 2.8. Suppose that F is a C 1 -function and x * is an R-regular solution of NCP(F). Then x * is a BD-regular solution of the system F (x) = 0.

A crucial property in the design of globally convergent algorithms, that has been extensively exploited in recent years, is that the corresponding merit function F (see ( 2)) is continuously differentiable [START_REF] Facchinei | A new merit function for nonlinear complementarity problems and a related algorithm[END_REF][START_REF] Geiger | On the resolution of monotone complementarity problems[END_REF][START_REF] Jiang | A new nonsmooth equations approach to nonlinear complementarity problems[END_REF]. Proposition 2.9. Suppose that F is a C 1 -function. Then the merit function F is continuously differentiable with gradient ∇ F (x) = H T F (x) for an arbitrary matrix H ∈ ∂ F (x).

In [START_REF] De Luca | A semismooth equation approach to the solution of nonlinear complementarity problems[END_REF] it has been shown how to calculate elements in the B-subdifferential of F , and this, in view of Proposition 2.9, also allows to easily evaluate the gradient of F .

We now consider analogous issues for P . To this end let x * be a solution of NCP(F) and let H ∈ ∂ B P (x * ) be arbitrarily chosen. Let H i. denote the ith row of this matrix. Furthermore, let α, β and γ denote the index sets defined in the previous subsection. From the very definition of the operator P and the B-subdifferential ∂ B P (x * ), we have the following:

(a) If i ∈ α, then H i. = ∇ F i (x * ) T . (b) If i ∈ γ , then H i. = e T i . (c) If i ∈ β, then H i. = ∇ F i (x * ) T or H i. = e T i .
Hence there is an index set δ such that ∅ ⊆ δ ⊆ β and

H = F (x * ) α∪δ,α∪δ F (x * ) α∪δ,γ ∪ δ 0 γ ∪ δ,α∪δ I γ ∪ δ,γ ∪ δ , (8)
where δ := β\δ denotes the complementary subset of δ in β. By using this expression of ∂ B P (x * ) and the very definition of b-regularity, it is then immediate to prove the following proposition.

Proposition 2.10. Suppose that F is a C 1 -function and x * is a b-regular solution of NCP(F). Then x * is a BD-regular solution of the system P (x) = 0.

It is interesting to note that P is BD-regular at a solution of the complementarity problem under an assumption which is weaker than that needed to establish an analogous result for F (even is we saw in the previous section that these two assumptions coincide in the case of P 0 -functions). One may wonder whether this gap is only due to a lack in the analysis of the properties of F ; the following example shows that the gap is actually intrinsic.

Example 2.1. Let n = 2 and F : IR 2 → IR 2 be defined by

F(x) := -x 1 + x 2 -x 2 .
Then it is easy to see that the complementarity problem NCP(F) has the unique solution

x * := (0, 0) which is b-regular (since F (x * ) is a nondegenerate matrix) but not R-regular.
Now consider the sequence {x k } defined by

x k := 1 k , 2 k T .
Then x k → x * and F is continuously differentiable at all x k with Jacobian

F (x k ) = 0 1 √ 2 -1 0 √ 2
for all k. Obviously, the sequence { F (x k )} converges to the singular matrix

0 1 √ 2 -1 0 √ 2 ∈ ∂ B F (x * ),
i.e., x * is not a BD-regular solution of the system F (x) = 0.

On the other hand it is very important to note that, contrary to F , the natural merit function associated to the operator P is not differentiable everywhere in general, making the design of globally convergent algorithms based on P fairly cumbersome.

Algorithm and convergence

In this section we propose a general algorithmic scheme which is similar to the algorithms considered, e.g., in [START_REF] De Luca | A semismooth equation approach to the solution of nonlinear complementarity problems[END_REF][START_REF] Jiang | A new nonsmooth equations approach to nonlinear complementarity problems[END_REF]. The only difference is that, instead of using a specific direction, we use a generic direction on which we impose certain assumptions. In the next section, we shall illustrate how directions satisfying these assumptions can be calculated. This allows us to analyze in a unified way different algorithms, corresponding to different directions.

In the sequel, in order to prove global convergence, we shall always assume that the search direction used in the algorithmic scheme satisfies Assumption 1 below. Assumptions 2 sup , 2 quad and 2 fin will be invoked only to establish convergence rates.

Assumption 1 (a) d k = 0 ⇒ ∇ F (x k ) = 0. (b) If x k → x * and d k → 0, then ∇ F (x * ) = 0.

Assumption 2 sup

Let x * be a solution of NCP(F) and assume that {x k } is a sequence converging to x * . Then

lim k→∞ x k + d k -x * x k -x * = 0.

Assumption 2 quad

Let x * be a solution of NCP(F) and assume that {x k } is a sequence converging to x * . Then

lim sup k→∞ x k + d k -x * x k -x * 2 < +∞.

Assumption 2 fin

Let x * be a solution of NCP(F). There exists a neighborhood of x * such that, if x k belongs to , then

x * = x k + d k .
Note that in our three different Assumptions 2, no connection is assumed between the way the sequences {x k } and {d k } are generated. The three Assumptions 2 sup , 2 quad and 2 fin say that the direction has to be locally superlinearly convergent, locally quadratically convergent or has a finite termination property, respectively.

In what follows, a general scheme of an algorithm for the solution of the complementarity problem using the merit function F is given.

General Line Search Algorithm 0. Choose x 0 ∈ IR n , ρ > 0, s > 1, β ∈ (0, 1/2), σ ∈ (0, 1). Set k := 0. 1. If x k satisfies a suitable termination criterion: STOP. 2. Calculate a direction d k satisfying Assumption 1 and set p k = d k . 3. If F (x k + p k ) ≤ σ F (x k ), (9)
set x k+1 := x k + p k , set k := k + 1 and go to Step 1. 4. If p k does not satisfy the following test

∇ F (x k ) T p k ≤ -ρ p k s , ( 10
)
set p k = -∇ F (x k ). Find the smallest i k ∈ {0, 1, 2, . . .} such that F x k + 2 -i k p k ≤ F (x k ) + β2 -i k ∇ F (x k ) T p k , ( 11
)
set x k+1 := x k + 2 -i k p k , set k := k + 1 and go to Step 1.
The algorithm itself, and the corresponding convergence analysis, are heavily based on the continuous differentiability of F described in Section 2.3.

In the sequel, in order to facilitate the statement of the convergence properties of the algorithm, we shall always assume that the termination criterion at Step 1 is never satisfied, so that an infinite sequence of points is generated Theorem 3.1. If Assumption 1 holds, then every limit point of a sequence {x k } generated by the General Line Search Algorithm is a stationary point of F .

Proof: The proof is by contradiction. Suppose, renumbering if necessary, that {x k } → x * and that ∇ F (x * ) = 0; then we can assume without loss of generality that the test ( 9) is never passed and that 0 < δ ≤ p k ≤ D [START_REF] Ferris | Engineering and economic applications of complementarity problems[END_REF] for suitable constants 0 < δ ≤ D. In fact, if the test ( 9) is satisfied infinitely many times, this would imply, recalling that at each step

F (x k+1 ) ≤ F (x k ), that { F (x k )} → 0, so that x * is a global minimum point of F and hence ∇ F (x * ) = 0. On the other hand if, for some subsequence K , { p k } K → 0, we have that ∇ F (x * ) = 0 by Assumption 1(b),
while { p k } cannot be unbounded because, taking into account that ∇ F (x k ) is bounded and s > 1, this would contradict [START_REF] Facchinei | A nonsmooth inexact Newton method for the solution of large-scale nonlinear complementarity problems[END_REF].

Then, since at each iteration [START_REF] Facchinei | A new merit function for nonlinear complementarity problems and a related algorithm[END_REF] holds and F is bounded from below, we have that { F (x k+1 ) -F (x k )} → 0 which implies, by the linesearch test,

2 -i k ∇ F (x k ) T d k → 0. ( 13
)
We want to show that 2 -i k is bounded away from 0. Suppose the contrary. Then, subsequencing if necessary, we have that {2 -i k } → 0 so that at each iteration the stepsize is reduced at least once and [START_REF] Facchinei | A new merit function for nonlinear complementarity problems and a related algorithm[END_REF] gives

F x k + 2 -(i k -1) p k -F (x k ) 2 -(i k -1) > β∇ F (x k ) T p k . ( 14
)
By [START_REF] Ferris | Engineering and economic applications of complementarity problems[END_REF] we can assume, subsequencing if necessary, that { p k } → p * = 0, so that, passing to the limit in ( 14), we get

∇ F (x * ) T p * ≥ β∇ F (x * ) T p *
and therefore

∇ F (x * ) T p * ≥ 0 ( 15
)
since β ∈ (0, 1). On the other hand, we also have, by [START_REF] Facchinei | A nonsmooth inexact Newton method for the solution of large-scale nonlinear complementarity problems[END_REF], that ∇ F (x * ) T p * ≤ -ρ p * s < 0, which contradicts (15); hence 2 -i k is bounded away from 0. But then [START_REF] Fischer | A special Newton-type optimization method[END_REF] and [START_REF] Facchinei | A nonsmooth inexact Newton method for the solution of large-scale nonlinear complementarity problems[END_REF]

imply that { p k } → 0 so that ∇ F (x * ) = 0 by Assumption 1(b). 2
Assumptions under which a stationary point of the function F is a solution of the complementarity problem are given in [START_REF] De Luca | A semismooth equation approach to the solution of nonlinear complementarity problems[END_REF] to which we refer the interested reader. Here we only note that if F is a P 0 -function then every stationary point of F is a solution of the complementarity problem. This condition is probably the weakest known condition guaranteeing the correspondence between stationary points of a merit function and solutions of the complementarity problem, and is one of the main motivations for using F instead of other merit functions in our General Line Search Algorithm.

In the sequel we examine the convergence rate of the algorithm. The proof uses the following two lemmas. The first one was shown in [START_REF] Facchinei | A new merit function for nonlinear complementarity problems and a related algorithm[END_REF] under a slightly stronger assumption. However, it is easy to see that the proof goes through with the weaker assumption stated here. The second lemma readily follows from [37, Lemma 3.1].

Lemma 3.3. Let G : IR n → IR n be semismooth and x * ∈ IR n be a BD-regular solution of the system G(x) = 0. Suppose that there are two sequences {x k } and {d k } such that

x k → x * and lim k→∞ x k + d k -x * x k -x * = 0.
Then

lim k→∞ G(x k + d k ) G(x k ) = 0.

Lemma 3.3.

There exist constants c 1 > 0 and c 2 > 0 such that

c 1 P (x) ≤ F (x) ≤ c 2 P (x)
holds for all x ∈ IR n , where P (x) := 1 2 P (x) 2 denotes the merit function belonging to the Minimum operator P .

Proof of Theorem 3.2: By Lemmas 3.3, 3.4, Proposition 2.10 and Assumption 2 sup , the b-regularity condition implies that eventually the test ( 9) is satisfied. So the assertion follows by the instructions of the General Line Search Algorithm and Assumption 2 sup .

2

The following two results can be easily proved along the same lines used for the proof of Theorem 3.2.

Theorem 3.5. Suppose that x * is a b-regular solution of NCP(F), that {x k } converges to x * and that Assumption 2 quad is satisfied in x * . Then, eventually p k = d k , x k+1 = x k + d k and {x k } converges Q-quadratically to x * . Theorem 3.6. Suppose that x * is a b-regular solution of NCP(F), that {x k } converges to x * and that Assumption 2 fin is satisfied in x * . Then there exists an index k 0 such that x k 0 = x * .

Search directions

By the results of the previous section we see that we can define algorithms which enjoy favourable properties if we are able to define search directions which satisfy Assumptions 1 and 2 sup (or 2 quad or 2 fin ). Obviously, from a theoretical point of view, we wish these assumptions to hold under conditions on the complementarity problem which are as weak as possible. In this section we consider several possible directions d k and compare their theoretical properties. A numerical comparison is the subject of Section 5.

The Fischer-Qi direction

The first direction we consider is basically the semismooth Newton direction for the system

F (x) = 0. d k = solution of H k d = -F (x k ), with H k ∈ ∂ B F (x k ), if system solvable (16) -∇ F (x k )
otherwise. Proof:

(a) The proof of this point is obvious if

d k = -∇ F (x k ). So consider the case in which d k is the solution of H k d = -F (x k ). If x k → x * (in particular, if x k = x * for every k) and if d k → 0,

The inexact Fischer-Qi direction

This is nothing else but the previous direction in which, however, the linear system is solved inexactly.

d k = solution of H k d = -F (x k ) + r k , with H k ∈ ∂ B F (x k ), if system solvable -∇ F (x k ) otherwise. ( 17
)
Theorem 4.2. Assume that the sequence of residuals

{r k } satisfies r k ≤ η k F (x k ) ,
where {η k } is a sequence of numbers such that, for every k, 0 ≤ η k ≤ η < 1. Then:

(a) The direction defined by [START_REF] Gomes-Ruggiero | Solving nonsmooth equations by means of quasi-Newton methods with globalization[END_REF]

The inexact LM Fischer-Qi direction

The direction considered in this subsection is an inexact perturbed Levenberg-Marquardt version of direction [START_REF] Geiger | On the resolution of monotone complementarity problems[END_REF]. It was already considered in [START_REF] Facchinei | A nonsmooth inexact Newton method for the solution of large-scale nonlinear complementarity problems[END_REF] and is defined as follows:

H T k H k + σ k I d = -H T k F (x k ) + r k , ( 18
)
where σ k is a nonnegative number and r k is the residual, which measures the inaccuracy with which the system is solved. Recall that

(H k ) T F (x k ) = ∇ F (x k
) and note that the linear system (18) always admits an exact (and hence inexact) solution since, for σ k > 0, the matrix on the left-hand side is positive definite, whereas for σ k = 0 we try to solve the normal equations of the system H k d = -F (x k ).

Theorem 4.3. Assume that the sequence of residuals {r

k } satisfies r k ≤ η k ∇ F (x k ) ,
where {η k } is a sequence of numbers such that, for every k, 0 ≤ η k ≤ η < 1. Assume also that the sequence {σ k } is bounded from above. Then: (a) The direction defined by [START_REF] Gomes-Ruggiero | Solving nonsmooth equations by means of quasi-Newton methods with globalization[END_REF]

satisfies Assumption 1. (b) If x * is an R-regular solution of NCP(F), r k = o( ∇ F (x k ) ) and {σ k } → 0, then Assumption 2 sup holds. Furthermore, if F is an LC 1 -function, r k = O( ∇ F (x k ) 2 ) and σ k = O( ∇ F (x k ) )
, then also Assumption 2 quad holds.

Proof:

(a) If x k → x * (in particular, if x k = x * for every k), we have, renumbering if necessary and taking into account the boundedness of the sequence {r k }, that r k → r * for some vector r * ∈ IR n . Hence, if d k → 0, we also have, using the boundedness of the Bsubdifferential on bounded sets and the assumed boundedness of the sequence {σ k }, that ∇ F (x * ) = r * . But, since r * ≤ η ∇ F (x * ) with η < 1 in view of our assumptions, this is only possible if ∇ F (x * ) = r * = 0. (b) If x * is an R-regular solution of NCP(F), then, by Proposition 2.8, x * is a BD-regular solution of the system F (x) = 0, which is semismooth by Proposition 2.7(a). Hence Assumption 2 sup holds because of Theorem 2.3(a) and (b). If, in addition, F is an LC 1 -function, then F is strongly semismooth by Proposition 2.7(b), and the assertion follows from Theorem 2.3(c). 2

The Pang-Qi direction

This search direction is basically the semismooth Newton direction for the system P (x) = 0. Point (c) in the theorem above also appears to be strongly related to the finite termination result in [START_REF] Facchinei | A new merit function for nonlinear complementarity problems and a related algorithm[END_REF].

d k = solution of H k d = -P (x k ), with H k ∈ ∂ B P (x k ), if system solvable -∇ F (x k ) otherwise. (19)
From a computational point of view, it may be interesting to observe that the system H k d = -P (x k ), used in [START_REF] Jiang | A new nonsmooth equations approach to nonlinear complementarity problems[END_REF], is very structured, and this can be exploited numerically. In fact, defining the index sets

α := α(x k ) := i | x k i > F i (x k ) , β := β(x k ) := i | x k i = F i (x k ) , γ := γ (x k ) := i | x k i < F i (x k )
and recalling the definition of P , we see, similarly to [START_REF] Dirkse | Crash techniques for large-scale complementairty problems[END_REF], that we can write the system

H k d = -P (x k ) as F (x k ) α∪δ,α∪δ F (x k ) α∪δ,γ ∪ δ 0 γ ∪ δ,α∪δ I γ ∪ δ,γ ∪ δ d α∪δ d γ ∪ δ = - F α∪δ (x k ) x k γ ∪ δ (20)
for some index set δ with ∅ ⊆ δ ⊆ β and δ := β\δ. This readily gives

d γ ∪ δ = -x k γ ∪ δ .
It is then obvious that it is not very sensible to apply a linear solver to the whole linear system in order to get a solution of [START_REF] Kanzow | Semismooth Newton-type methods for the solution of nonlinear complementarity problems[END_REF]. Is is more convenient to set

d k γ ∪ δ = -x k γ ∪
δ , substitute this value in the first set of equations in [START_REF] Kanzow | Semismooth Newton-type methods for the solution of nonlinear complementarity problems[END_REF] and solve the resulting reduced linear system, which reads

F α∪δ,α∪δ (x k )d α∪δ = -F α∪δ (x k ) + F (x k ) α∪δ,γ ∪ δ x k γ ∪ δ .
(21)

The inexact Pang-Qi direction

This is nothing else but the previous direction in which, however, the linear system is solved approximately.

d k = solution of H k d = -P (x k ) + r k , with H k ∈ ∂ B P (x k ), if system solvable -∇ F (x k ) otherwise. ( 22
)
The proof of the following theorem, which we omit, can be carried out exactly along the same lines used in Theorem 4. Obviously, also in this case, similarly to what has been seen in the previous section, it is possible to approximately solve just the reduced system (21).

The inexact LM Pang-Qi direction

The direction considered in this subsection is an inexact perturbed Levenberg-Marquardt version of direction [START_REF] Jiang | A new nonsmooth equations approach to nonlinear complementarity problems[END_REF].

H T k H k + σ k I d = -H T k P (x k ) + r k , ( 23
)
where H k ∈ ∂ B P (x k ), σ k is a nonnegative number and r k is the residual which measures the inaccuracy with which the system is solved.

Theorem 4.6. Assume that the sequence of residuals {r k } satisfies r k ≤ η k H T k P (x k ) , where {η k } is a sequence of numbers such that, for every k, 0 ≤ η k ≤ η < 1, and assume that the sequence {σ k } is bounded from above. Assume also that for every k the matrix F α∪δ,α∪δ (x k ) is nonsingular (where α and δ denote the index sets used in Subsection 4.4). Then: (a) The direction defined by [START_REF] Gomes-Ruggiero | Solving nonsmooth equations by means of quasi-Newton methods with globalization[END_REF]

satisfies Assumption 1. (b) If x * is a b-regular solution of NCP(F), r k = o( H T k P (x k ) ) and {σ k } → 0, then Assumption 2 sup holds. Furthermore, if F is an LC 1 -function, r k = O( H T k P (x k ) 2 ) and σ k = O( H T k P (x k )})
, then also Assumption 2 quad holds.

Proof:

(a) If x k → x * (in particular, if x k = x * for every k), we have, renumbering if necessary and taking into account the boundedness of the sequence {r k }, that r k → r * for some vector r * ∈ IR n . Hence, if d k → 0, we also have, using the boundedness of the B-subdifferential on bounded sets, the boundedness of the sequence {σ k } and the fact that the B-subdifferential is a closed mapping, that H T * P (x * ) = r * for some matrix H * ∈ ∂ B P (x * ). In view of our assumptions, however, we also have r * ≤ η H T * P (x * ) with η < 1. But this is only possible if H T * P (x * ) = r * = 0. In turn, recalling [START_REF] Kanzow | Semismooth Newton-type methods for the solution of nonlinear complementarity problems[END_REF] and the nonsingularity assumption, this implies P (x * ) = 0. Hence x * is a solution of NCP(F). This implies Also in this case, as seen in the previous two cases, things can be arranged so that only a reduced system has to be solved approximately at each iteration. Note also that, in order to satisfy Assumption 1, we had to impose a nonsingularity assumption in Theorem 4.6.

∇ F (x * ) = 0. (b) If x * is a b-regular solution of NCP(F)
This assumption, which is satisfied in a sufficiently small neighbourhood of a b-regular solution of NCP(F), was not necessary for any of the other search directions discussed in this section.

Comparison of the theoretical characteristics of the directions

From a theoretical point of view, the Pang-Qi directions (exact, inexact and Levenberg-Marquardt) need weaker assumptions than the corresponding Fischer-Qi directions in order to satisfy Assumption 2 sup or Assumption 2 quad , namely b-regularity instead of R-regularity.

However it should also be added that, as shown in Section 2.3, these two conditions are equivalent for the class of P 0 -functions. Another advantage of the Pang-Qi directions is that they only need the (possibly inexact) solution of a linear system of a dimension that, asymptotically, is equal to the number of positive variables at the solution. This contrasts favourably with the Fischer-Qi directions that always need the solution of an n-dimensional square system. Finally, it should also be noted that there is only one direction for which Assumption 2 fin (i.e., finite termination for linear complementarity problems) can be proved: the exact Pang-Qi direction. The Fischer-Qi directions seem to have only one, but important, advantage over the corresponding Pang-Qi directions: We saw that, in the General Line Search Algorithm, we first try to use a Newton-like direction, but, if neither (9) nor [START_REF] Facchinei | A nonsmooth inexact Newton method for the solution of large-scale nonlinear complementarity problems[END_REF] are satisfied, we employ -∇ F (x k ) as search direction. Obviously, it does not appear to be very desirable to have to resort to the gradient often. From this point of view we may note that it is easy to see that the Fischer-Qi directions are always directions of descent for the merit function F , see, e.g., [START_REF] De Luca | A semismooth equation approach to the solution of nonlinear complementarity problems[END_REF][START_REF] Facchinei | A new merit function for nonlinear complementarity problems and a related algorithm[END_REF], while we cannot expect a similar property to hold for the Pang-Qi directions since these directions are calculated on the basis of a reformulation of the complementarity problem based on P and not on F . Thus, a priori, we should expect that, if we use a Pang-Qi direction in the General Line Search Algorithm, the Newton direction will be discarded more often than when we use a Fischer-Qi direction.

We conclude this section by noting that some of the algorithms obtained by using some of the directions introduced in this section in the General Line Search Algorithm have already been presented in the literature. In particular, if we use the exact Fischer-Qi direction we have the algorithm considered in [START_REF] De Luca | A semismooth equation approach to the solution of nonlinear complementarity problems[END_REF], whereas if we use the inexact LM Fischer-Qi direction we obtain the algorithm considered in [START_REF] Facchinei | A nonsmooth inexact Newton method for the solution of large-scale nonlinear complementarity problems[END_REF].

Numerical results

Description of the test problems and implementation details

We tested our algorithm on test problems generated by using a technique first suggested in [START_REF] Gomes-Ruggiero | Solving nonsmooth equations by means of quasi-Newton methods with globalization[END_REF] and subsequently used, e.g., in [START_REF] Facchinei | A nonsmooth inexact Newton method for the solution of large-scale nonlinear complementarity problems[END_REF]. Let g(x) = 0 be a (large-scale) differentiable system of nonlinear equations and let x * ∈ IR n be defined by x * = (1, 0, 1, 0, . . .) T . For all i = 1, . . . , n set

F i (x) = g i (x) -g i (x * ) if i odd or i > r, g i (x) -g i (x * ) + 1 otherwise,
where r ≥ 0 is a given integer. Then x * is a solution of the nonlinear complementarity problem NCP(F) (but not necessarily its unique solution). Note that, if r is smaller than n -1, the problem is degenerate at x * . As done in [START_REF] Gomes-Ruggiero | Solving nonsmooth equations by means of quasi-Newton methods with globalization[END_REF], we used the collection of 17 large-scale problems from Lukšan [START_REF] Lukšan | Inexact trust region method for large sparse systems of nonlinear equations[END_REF] and the starting points indicated there. However, we did not consider problems 3 and 6, since these were never solved by any of the algorithms considered; actually we are not aware of any algorithm capable of solving these problems (see also [START_REF] Gomes-Ruggiero | Solving nonsmooth equations by means of quasi-Newton methods with globalization[END_REF]).

We considered problems with dimensions n = 100, n = 1000 and n = 10000. For each dimension we considered two cases: r = n/2 and r = n. The former case corresponds to a degenerate problem with n/4 degenerate components, while the latter case corresponds to a nondegenerate problem. Finally, besides the starting points x 0 suggested in [START_REF] Lukšan | Inexact trust region method for large sparse systems of nonlinear equations[END_REF], we also considered, for each problem, an additional starting point defined by

x0 i = 10x 0 i if x 0 i = 0, 10 otherwise.
Therefore each algorithm was tested on 90 different problems and on each of these problems two different starting points were used so that we have 180 runs for each algorithm.

In the previous section, we introduced 6 directions which, coupled with the General Line Search Algorithm, give 6 different algorithms. However, in this section we only report the results for four directions: the exact Fischer-Qi (EXFQ) direction, the exact Pang-Qi (EXPQ) direction, the Levenberg-Marquardt inexact Fischer-Qi (LMFQ) direction and the Levenberg-Marquardt inexact Pang-Qi (LMPQ) direction. We do not report results on the inexact Fischer-Qi and inexact Pang-Qi directions. In fact, in order to implement the corresponding algorithms we need an iterative solver for unsymmetric linear systems. Although many of them are available, our experiments seem to indicate that these solvers (at least if they are not combined with a sophisticated preconditioner) are not very reliable in practice (at least on our problems) leading either to extremely long running times or to failure due to the incapability of finding a sufficiently good search direction. We were genuinely surprised by these results. We think that the problem is that these solvers are usually thought for and tested on linear systems arising from the discretization of PDE problems, and they are not very reliable when used on other problems. For the other 4 methods, instead, we had no difficulty in selecting an appropriate solver. In particular, we used HARWELL routine MA50 to calculate the EXFQ and EXPQ directions and the conjugate gradient method to calculate the LMFQ and LMPQ directions.

The conjugate gradient algorithm was stopped when the norm of the residual is smaller than (0.1/(k + 1)) ∇ F (x k ) in the case of the LMFQ algorithm and when the residual is smaller than

0.1 k + 1 H k α∪δ,α∪δ T F α∪δ (x k ) -F (x k ) α∪δ,γ ∪ δ x k γ ∪ δ
in the case of the LMPQ algorithm. We also set a limit of 200 iterations to the conjugate gradient phase. Finally we turn to the choice of σ k . We first note that the classical sophisticated choices indicated, e.g., in [START_REF] Dennis | Numerical Methods for Unconstrained Optimization and Nonlinear Equations[END_REF], are not suitable for large-scale problems. We therefore used a very simple, and yet seemingly effective strategy, already employed in [START_REF] Facchinei | A nonsmooth inexact Newton method for the solution of large-scale nonlinear complementarity problems[END_REF]. If in the previous iteration the quotient ∇ F (x k-1 ) / d k-1 is greater than 250 and the norm of the natural residual min{x k , F(x k )} is greater than k(0.1 √ n) then we set σ k = 1, otherwise we set σ k = 0. The first test is a rough indicator that "something is going wrong" while the second test makes the possibility of the perturbation (i.e., σ k > 0) more and more unlikely the more the process progresses, so that the final fast convergence rate is preserved.

We coded the algorithm in Fortran 77 and run it on an IBM RISC 6000/375 machine in double precision arithmetic.

The main stopping criterion is min{x k , F(x k )} ≤ 10 -5 √ n or ∇ F (x k ) ≤ 10 -5 √ n, but we also stopped the algorithm after 100 iterations if the former stopping criterion was not met. We used a single set of parameters for all the runs of the problems, more precisely we set: ρ = 10 -8 , s = 2.1, β = 10 -4 and σ = 0.9.

Finally, we also considered the option of initiating the computation by performing at most 10 iterations of a projected gradient method to minimize F over the nonnegative orthant [START_REF] Bertsekas | Nonlinear Programming[END_REF]. The rationale behind this option is to try to move in a more "promising" zone by performing a few steps of a robust and not expensive method. Similar ideas were successfully used, for example, in [START_REF] Dirkse | Crash techniques for large-scale complementairty problems[END_REF], where, however, the situation is slightly more complicated since the merit function is not differentiable. It should be remarked that in this gradient phase, we introduce the constraints x ≥ 0. This seems reasonable, since, on the one hand, we know that the solution of the problem is in the nonnegative orthant and, on the other hand, the projected gradient method can effectively handle the bound constraints. Since it is known that the projected gradient method usually behaves well in the first iterations and then becomes slow, we impose a limit of ten iterations to this phase. However, we stopped earlier if one of the following criteria is satisfied:

-( F (x k-1 ) -F (x k ))/ F (x k ) ≤ 0.05;
-the zero variables are the same in x k-1 and x k and

( F (x k-1 ) -F (x k ))/ F (x k ) ≤ 0.1; -F (x k ) ≤ 10 -5 √ n.

Analysis of the results

Our first aim is to determine whether the initial steps of projected gradient are beneficial.

In figure 1 we report the number of failures of each algorithm without using the projected gradient phase (basic procedure) and with the projected gradient phase (modified procedure). It is apparent that the projected gradient phase greatly enhances the robustness of all the algorithms. The overall failures pass from 140 to 72, with a reduction of 50%. We also see that the algorithms which seem to benefit more by the gradient steps are the exact ones. This probably indicates that a pure Newton direction is less robust than a mixture of Newton and gradient directions. In fact, the Levenberg-Marquardt directions may be seen as a way of combining the Newton and gradient direction. Nevertheless, even this seems to be not enough and a few initial gradient steps improve the performance also of the Levenberg-Marquardt algorithms. It may also be noted that if the gradient projection strategy is used, there seems to be no significant difference in the robustness of the four algorithms, thus suggesting that the gradient projection phase is actually effective in moving the starting point in a good region, from which all the four directions considered behave well.

In figure 2 we report more in detail the number of failures for each algorithm and for each set of test problems. Some observations can be made. The degeneracy of the solution, that corresponds to the nonsmoothness of F and P at the solution, seems to have no impact on the robustness of the algorithms. Furthermore, as could easily be expected, the faraway starting points are consistently more difficult than the original ones. What is more surprising is the relatively high number of failures for the exact algorithms (without projected gradient phase), even for the standard starting points. Again, this seems to suggest that the pure Newton directions are not very robust, and that taking into account the gradient direction, as is done by the Levenberg-Marquardt methods, increases the stability of the algorithms.

We think that the previous results clearly show that the gradient phase is overall very useful in increasing the robustness of all algorithms considered, even if it appears to be more useful in the case of the exact algorithms. The next step is then to determine how costly the projected gradient phase is and which algorithm is more efficient from the computational point of view. In order to compare the various algorithms in a homogeneous way, we only considered the subset of test problems for which all algorithms (i.e., the four basic algorithms and the ones including the gradient phase) did not fail and converged to the same solution. This results in a total of 55 test problems: 22 for n = 100, 16 for n = 1000 and 17 for n = 10000.

In figures 3-5 we report the cumulative times needed by the algorithms to solve the problems of dimension n = 100, n = 1000 and n = 10000, respectively. We reported the times according to the dimensions because the times needed to solve problems, let us say, of dimension 100 are not comparable to those needed to solve those of dimension 10000. Let us first compare the version using exact linear solvers (EXPQ and EXFQ). We see that the modified procedure, i.e., the one using the preliminary steps of projected gradient, never requires more time than the basic procedure. This fact, along with the improvement in the number of failures, seems to indicate that the projected gradient phase is very beneficial on these two algorithms. Another observation is that the algorithms using the Pang-Qi direction are consistently faster than their counterparts using the Fischer-Qi direction. This behavior is due to the fact that to calculate the Pang-Qi direction it is sufficient to solve a linear system of reduced dimension, as already observed in Section 4.4.

If we pass to examine the inexact Levenberg-Marquardt algorithms (LMPQ and LMFQ), we see that some differences are present. In particular, the time needed by the modified versions can be substantially higher than the one needed by the corresponding basic algorithm. In particular, this is apparent in the case n = 1000 for both algorithms and in the case n = 10000 for the LMPQ one. A closer look at the runs showed that this anomalous behavior is almost completely due to the behavior on three instances of problem 15. When solving these test problems, the projected gradient phase passes an initial point to the Newton phase that gives rise to a Levenberg-Marquardt system which is extremely difficult to solve by the conjugate gradient method. The increase of times reported in figures 4 and 5 is therefore almost entirely due to the very high number of conjugate gradient inner iterations needed to approximately solve the first systems of the Newton phase when solving these instances of problem 15.

In figure 6 and in figure 7 we therefore report the same data of figures 4 and 5 without considering, however, the three instances of problem 15. If we now look at the figures 3, 6 and 7, we see that, similarly to what observed for the EXPQ and EXFQ algorithms, the LMPQ algorithms are faster than the LMFQ ones. If we compare the basic and the modified inexact Levenberg-Marquardt versions, we see that the modified versions usually need the same time as the basic versions (except that for the LMFQ when n = 10000). This is remarkably different from what we observed for the exact algorithms, where a substantial decrease of the time takes place.

To understand better this behavior, we report in figure 8 the number of iterations of the various algorithms. Note that the number of iterations of the projected gradient is, obviously, the same for every algorithm, and that the effect of the gradient phase is to reduce the number of subsequent Newton-type iterations (compare the black columns with the dark grey ones). The percentual decrease in the number of Newton iterations, when passing from the basic procedure to the modified one, is roughly similar for all the algorithms. It is then natural that the weight of the gradient phase will be higher for the Levenberg-Marquardt algorithms than for the exact ones since, at least on the problems considered here, the work needed to solve exactly one linear system is higher than the cost to approximately solve the same system by the conjugate gradient method (recall that one system is solved, possibly approximately, at each iteration).

With regard to figure 8, it is also interesting to note that the number of Newton iterations for the four algorithms (both in the basic or in the modified procedure) is very similar for all the algorithms. Thus we may conclude that the more marked differences in the times reported in the previous tables are mainly due to the different cost of the search direction computations and, possibly, to the number of function evaluations. In fact we recall that the number of Jacobian evaluations is equal to the number of total iterations, and therefore roughly the same for all algorithms. In figure 9 we therefore report the number of function evaluations. There are two interesting facts to be observed. First the high number function evaluations needed by the projected gradient phase. In general the projected gradient method is considered to be a cheap method, however, this is true only if the cost of function evaluations is "low". In our case low should be read as (much) cheaper than the cost of (approximately) solving a linear system. In our test problem set the function evaluations are actually very cheap; however, should the algorithms be applied to problems with very high function evaluation costs, the times of the modified procedures could become much higher than the cost of the basic procedure, contrary to what is reported in figures 3-7.

A second point to be remarked is that, both in the case of the basic and of the modified versions, the number of function evaluations needed in the Newton phase is sensibly lower for the algorithms based on the Fischer-Qi directions. This corresponds to the fact that much less backtrackings occur when using a Fischer-Qi direction than when using a Pang-Qi direction. This confirms the observation made in Section 4.7 that the Fischer-Qi direction is much more directly related to the merit function F while the Pang-Qi directions do not even need to be descent directions for F .

Again, since in our test set the cost of function evaluations is low, this behavior does not affect much the times reported in figures 3-7. However, we should expect that the EXFQ and LMFQ methods become more and more competitive as they are applied to problems with higher and higher function evaluation costs.

Summary of the results

The previous observations can be summarized as follows. The projected gradient phase certainly the robustness of all the four algorithms. The additional cost of this phase is low on the test problems considered here, but could become high for those problems with a high function evaluation cost. If the projected gradient phase is used, the reliability of the four algorithms is similar, and they roughly need the same number of iterations. If, however, we consider the times needed by the various algorithms, the behavior may differ widely. If the cost of the function evaluations is low, the Pang-Qi algorithms are probably more convenient; however, if we want to solve problems with functions extremely expensive to evaluate, the Fischer-Qi algorithms can be more attractive.

Another aspect that should be kept in mind is the difficulty in solving the linear systems. All the problems we considered gave rise to sparse, "easy" systems; and this is typical of many applications. If the linear systems are "difficult" and not huge, the exact versions could be preferable; on the other hand, if extremely large instances have to be solved, the inexact algorithms are probably the only available options.

Final remarks

We have presented and studied a general line search scheme which allows us to analyze in a unified framework several semismooth algorithms for the solution of nonlinear complementarity problems. Several search directions have been studied in detail, both from a theoretical and a numerical point of view. However, we remark that it is easily possible to define analyze in a similar way other algorithms. Below, we hint at some of these further possibilities.

We first recall that a mapping ϕ : IR We denote the corresponding operators by KK and CCK , respectively. Many other NCPfunctions can be contructed, see Mangasarian [START_REF] Mangasarian | Equivalence of the complementarity problem to a system of nonlinear equations[END_REF] as well as Sun and Qi [START_REF] Sun | On NCP-functions[END_REF].

The function ϕ KK is interesting because for λ = 2 it reduces to the Fischer function, and in the limiting case λ = 0 it becomes a multiple of the Minimum function. On the other hand, very strong theoretical and numerical results have been reported for ϕ CCK in [START_REF] Chen | A penalized Fischer-Burmeister function: Theoretical investigation numerical results[END_REF]. Now, similar to the (exact, inexact and LM inexact) Fischer-Qi directions, we can define (exact, inexact and LM inexact) search directions based on KK and CCK . Since it is known [START_REF] Kanzow | A new class of semismooth Newton-type methods for nonlinear complementarity problems[END_REF][START_REF] Chen | A penalized Fischer-Burmeister function: Theoretical investigation numerical results[END_REF] that KK and CCK are (strongly) semismooth if F is a C 1 -function (LC 1 -function) and that an R-regular solution of NCP(F) is a BD-regular solution for the nonlinear systems of equations KK (x) = 0 and CCK (x) = 0, it is not difficult to see that Theorems 4.1, 4.2 and 4.3 also hold for the corresponding search directions based on KK and CCK .

To conclude, we also mention that the results of this paper can be generalized to variational inequalities with box constraints. In fact, the method introduced very recently in [START_REF] Kanzow | Solving box constrained variational inequality problems by using the natural residual with D-gap function globalization[END_REF] can be viewed as an extension of one of the methods discussed in this paper.

Theorem 2 . 5 .

 25 Let x * ∈ IR n be a solution of NCP(F). Then the following statements are equivalent: (a) x * is an R-regular solution of NCP(F). (b) x * is a b-regular solution of NCP(F) and the Schur-complement (7) is a P-matrix. (c) x * is a b-regular solution of NCP(F) and the Schur-complement (7) is a P 0 -matrix. Proof: For simplicity, let us write M := F (x * ).

Proposition 2. 7 .

 7 (a) If F is a C 1 -function, then F and P are semismooth. (b) If F is an LC 1 -function, then F and P are strongly semismooth.

Theorem 3 . 2 .

 32 Suppose that x * is a b-regular solution of NCP(F), that {x k } converges to x * and that Assumption 2 sup is satisfied in x * . Then, eventually p k = d k , x k+1 = x k + d k and {x k } converges Q-superlinearly to x * .

Theorem 4. 1 .

 1 (a) The direction defined by[START_REF] Geiger | On the resolution of monotone complementarity problems[END_REF] satisfies Assumption 1. (b) If x * is an R-regular solution of NCP(F), then Assumption 2 sup holds. Furthermore, ifF is an LC 1 -function, then also Assumption 2 quad holds.

  we have, taking into account the boundedness of the B-subdifferential on bounded sets, that F (x * ) = 0. But then x * is a global minimum of F so that ∇ F (x * ) = 0. (b) If x * is an R-regular solution of NCP(F) then, by Proposition 2.8, x * is a BD-regular solution of the system F (x) = 0, which is semismooth by Proposition 2.7(a). Hence Assumption 2 sup holds because of Theorem 2.1. If, in addition, F is an LC 1 -function, then F is strongly semismooth by Proposition 2.7(b), and, once again, the assertion follows from Theorem 2.1. 2

Theorem 4. 4 .

 4 (a) The direction defined by[START_REF] Jiang | A new nonsmooth equations approach to nonlinear complementarity problems[END_REF] satisfies Assumption 1. (b) If x * is a b-regular solution of NCP(F), then Assumption 2 sup holds. Furthermore, if F is an LC 1 -function, then also Assumption 2 quad holds. (c) If F is affine and x * is a b-regular solution, then Assumption 2 fin holds.Proof:(a) The proof of this point is obvious ifd k = -∇ F (x k ). So consider the case in which d k is the solution of H k d = -P (x k ). If x k → x * (in particular, if x k = x * for every k) and if d k → 0,we have, taking into account the boundedness of the B-subdifferential on bounded sets, that P (x * ) = 0. But then x * is a solution of NCP(F) so that ∇ F (x * ) = 0. (b) If x * is a b-regular solution of NCP(F) then, by Proposition 2.10, x * is a BD-regular solution of the system P (x) = 0, which is semismooth by Proposition 2.7(a). Hence Assumption 2 sup holds because of Theorem 2.1(a) and (b). If, in addition, F is an LC 1 -function, then P is strongly semismooth by Proposition 2.7(b), and the assertion follows from Theorem 2.1(c). (c) This follows immediately from [15, Theorem 4], Proposition 2.10 and the upper semicontinuity of the B-subdifferential. 2

Figure 1 .

 1 Figure 1. Failures.

Figure 2 .

 2 Figure 2. Failures in detail.

Figure 3 .

 3 Figure 3. Times in 10 -2 × seconds, n = 100.

Figure 4 .

 4 Figure 4. Times in 10 -2 × seconds, n = 1000.

Figure 5 .

 5 Figure 5. Times in 10 -2 × seconds, n = 10000.

Figure 6 .

 6 Figure 6. Times in 10 -2 × seconds, n = 1000, without Problems 15.

Figure 7 .

 7 Figure 7. Times in 10 -2 × seconds, n = 10000, without Problems 15.

Figure 8 .

 8 Figure 8. Number of iterations.

Figure 9 .

 9 Figure 9. Number of function evaluations.

2 →

 2 IR is called an NCP-function if ϕ(a, b) = 0 ⇐⇒ a ≥ 0, b ≥ 0, ab = 0, cf. (1). If ϕ is any NCP-function, then NCP(F) is equivalent to (x) = 0,where : IR n → IR n is defined by(x) := ⎛ ⎜ ⎜ ⎝ ϕ(x 1 , F 1 (x)) . . . ϕ(x n , F n (x))Fischer function ϕ F and the Minimum function ϕ P are NCP-functions. Other examples include, e.g., ϕ KK (a, b) := (ab) 2 + λabab, λ ∈ (0, 4) fixed (see Kanzow and Kleinmichel [22]), and ϕ CCK (a, b) := λϕ F (a, b) -(1 -λ)a + b + , λ ∈ (0, 1) fixed (see Chen, Chen and Kanzow [2]), where a + := max{0, a} and, similarly, b + := max{0, b}.

  If the sequence {x k } converges to the solution x * and if G is strongly semismooth, then the rate of convergence is Q-quadratic if and only if r k = O( G(x k ) 2 ).

	* .
	(b) If the sequence {x k } converges to the solution x * , then the rate of convergence is Q-
	superlinear if and only if r k = o( G(x k ) ).
	(c) Paralleling well-known results in the smooth case, we can also consider a Levenberg-
	Marquardt type version of the Generalized Newton Method.
	Generalized Inexact Levenberg-Marquardt Method
	0. Choose x 0 ∈ IR n and set k := 0.
	1. Select an element H k ∈ ∂ B G(x k ) and a nonnegative σ k . Compute a solution d k ∈ IR n
	of the equation

  satisfies Assumption 1. (b) If x * is an R-regular solution of NCP(F) and r k = o( F (x k ) ), then Assumption 2 sup holds. Furthermore, if F is an LC 1 -function and r k = O( F (x k ) 2 ), then also Assumption 2 quad holds. The proof of this point is obvious ifd k = -∇ F (x k ). So consider the case in which d k is the solution of H k d = -F (x k ) + r k for a suitable residual vector r k ∈ IR n . If x k → x * (in particular, if x k =x * for every k), we have, renumbering if necessary and taking into account the boundedness of the sequence {r k }, that r k → r * for some vector r * ∈ IR n . Hence, if d k → 0, we also have, taking into account the boundedness of the Bsubdifferential on bounded sets, that F (x * ) = r * . In view of our assumptions, however, we have r * ≤ η F (x * ) . Since η < 1, this is only possible if F (x * ) = r * = 0. In turn, this implies that x * is a global minimum of F so that ∇ F (x * ) = 0. (b) If x * is an R-regular solution of NCP(F) then, by Proposition 2.8, x * is a BD-regular solution of the system

	Proof:
	(a)

F (x) = 0, which is semismooth by Proposition 2.7(a). Hence Assumption 2 sup holds because of Theorem 2.2(a) and (b). If, in addition, F is an LC 1 -function, then F is strongly semismooth by Proposition 2.7(b), and the assertion follows from Theorem 2.2(c).

2

  If x * is a b-regular solution of NCP(F) and r k = o( P (x k ) ), then Assumption 2 sup holds. Furthermore, if F is an LC 1 -function and r k = O( P (x k ) 2 ), then also Assumption 2 quad holds.

3, using Proposition 2.10 instead of Proposition 2.8. Theorem 4.5. Assume that the sequence of residuals {r k } satisfies r k ≤ η k P (x k ) , where {η k } is a sequence of numbers such that, for every k, 0 ≤ η k ≤ η < 1. Then: (a) The direction defined by (17) satisfies Assumption 1. (b)

  LC1 -function, then P is strongly semismooth by Proposition 2.7(b), and the assertion follows from Theorem 2.3(c).2

then, by Proposition 2.10, x * is a BD-regular solution of the system P (x) = 0, which is semismooth by Proposition 2.7(a). Hence Assumption 2 sup holds because of Theorem 2.3(a) and (b). If, in addition, F is an


Ibtissem Ben Aïcha

Youssef Mejri

SIMULTANEOUS DETERMINATION OF THE MAGNETIC FIELD AND THE ELECTRIC POTENTIAL IN THE SCHRÖDINGER EQUATION BY A FINITE NUMBER OF BOUNDARY OBSERVATIONS

Keywords: Inverse problem, magnetic Schrödinger equation, stability estimate, finite number of partial Neuman measurements, Carleman estimate

We study the inverse problem of determining the magnetic field and the electric potential appearing in the magnetic Schrödinger equation, from the knowledge of a finite number of lateral observations of the solution. We prove a Lipschitz stability estimate for both coefficients simultaneasly by choosing the "initial" conditions suitably.

INTRODUCTION

1.1. Statement of the problem. In this paper, we study the following inverse problem: Given T ¡ 0 and a bounded domain Ω R n , n ¥ 1, with smooth boundary Γ, we want to determine simultaneously the divergence free magnetic potential of the form apx, tq : χptqapxq and the electric potential qpx, tq : βptqqpxq appearing in the following equation (1.1) 6 8 7 p¡if t H a ptq qpx, tqqupx, tq 0, in Q Ω ¢ p0, T q, upx, T 2 q u 0 pxq in Ω, upx, tq 0, on Σ Γ ¢ p0, T q where H a ptq : pi∇ χptqaq 2 denotes the time-dependent Hamiltonian, associated to the magnetic potential vector χptqapxq. Here a pa 1 , ..., a n q A : H 1 pΩq n ta L V pΩ, R n q, ∇ ¤ a 0u, and q L V pΩq are unknown real valued functions. Moreover, the functions β, χ C 3 p0, T ; Rq are assumed to be known functions satisfying (1.2) χp T 2 q βp T 2 q 0, χ I p T 2 q $ 0, β I p T 2 q $ 0.

We denote by Γ an open subset of Γ satisfying an appropriate geometrical condition we shall precise later.

The inverse problem we investigate in this paper, is to know whether the knowledge of a finite number of Neumann measurements f ν u |Γ ¢p0,Tq : p∇u ¤ νq |Γ ¢p0,Tq uniquely determines apxq and qpxq simultaneously. Here νpxq denotes the unit outward normal to Γ at x.

Bibliographical notes.

There is a wide mathematical literature dealing with uniqueness and stability in inverse coefficient problems related to partial differential equations. In recent years, significant progress have been made in the recovery of magnetic potentials in magnetic Schrödinger equations from the Dirichletto Neuman map. As it was noted in [START_REF] Eskin | Inverse problems for the Schödinger equations with time-dependent electromagnetic potentials and the Aharonov-Bohn effect[END_REF], the Dirichlet-to Neuman map is invariant under a gauge transfromation of the magnetic potential. Namely, given φ C 1 pΩq such that φ |Γ 0, we have Λ a ∇φ Λ a . So, due to this obstruction to uniqueness, the best one could recover from the Dirichlet-to-Neumann map is the magnetic field da, where da is the exterior derivative of a interpreted as the one form °n j1 a j dx j .

The inverse problem of determining the magnetic field da and the electric potential q from boundary observations was first considered by Sun [START_REF] Sun | An inverse boundary value problem for the Schrödinger operator with vector potentials[END_REF], in the case n ¥ 3. He showed that da and q can be uniquely determined when a W 2,V , q L V and da is small in the L V norm. In [START_REF] Nakamura | Global identifiability for an inverse problem[END_REF], a uniqueness result was proved for C V magnetic potentials. In [START_REF] Eskin | Inverse problems for the Schödinger equations with time-dependent electromagnetic potentials and the Aharonov-Bohn effect[END_REF], Eskin proved a uniqueness result for an inverse problem for the Schrödinger equation with electromagnetic potentials, modulo a gauge transform of the recovery of the potentials from the Dirichlet-to Neuman map. In a recent work, Bellassoued and Choulli [START_REF] Bellassoued | Stability estimate for an inverse problem for the magnetic Schödinger equation from the Dirichlet-to-Neumann map[END_REF] considered the problem of recovering the magnetic potential da from the Dirichlet-to-Neumann map associated to the Schrödinger equation and proved in dimension n ¥ 2 a stability estimate of Hölder type.

In the Riemannian case, Bellassoued [2] proved recently a Hölder-type stability estimate in the recovery of the magnetic field da and the electric potential q from the knowledge of the Dirichlet-to-Neumann map associated to the Schrödinger equation with zero initial data. In the absence of the magnetic potential a, the problem of recovering the electric potential q on a compact Riemannian manifold was solved by Bellassoued and Dos Santos Ferreira [START_REF] Bellassoued | Stable determination of coefficients in the dynamical anisotropic Schrödinger equation from the Dirichlet-to-Neumann map[END_REF].

The problem of determining the magnetic field by a local Dirichlet-to Neuman map was solved by Dos Santos Ferreira, Kenig, Sjöstrand and Uhlmann [START_REF] Santos Ferreira | Determining a magnetic Schrödinger operator from partial Cauchy data[END_REF]. In [START_REF] Tzou | Stability estimates for coefficients of magnetic Shrödinger equation from full and partial boundary measurements[END_REF], Tzou showed that the magnetic field depends stably on the Dirichlet-to Neuman map measured on any subboundary which is larger than half the boundary. In [START_REF] Joud | Stability estimate for an inverse problem for the Schrödinger equation in a magnetic field from partial boundary meausurements[END_REF], Benjoud studied the inverse problem of recovering the magnetic field da and the electric potential q from the knowledge of the Dirichlet-to-Neumann map. Assuming that the potentials are known in a neighborhood of the boundary, she proved a stability estimate with respect to arbitrary partial boundary observations. The key ingredient in the proof of all the above mentioned papers, is the construction of geometric optics solutions.

To our knowledge, there is a few results on the recovery of coefficients appearing in a Schrödinger equation, from a finite number of boundary measurements. By a method based essentially on an appropriate Carleman estimate, Baudouin and Puel [START_REF] Baudouin | Uniqueness and stability in an inverse problem for the Schrödinger equation[END_REF] showed that the electric potential in the Shrödinger equation can be stably recovered from a single boundary measurement. In [START_REF] Cristofol | Stability estimate in an inverse problem for non autonomous magnetic Schrödinger equations[END_REF], Cristofol and Soccorsi proved a Lipschitz stability in recovering the magnetic field in the Schrödinger equation from a finite number of observations, measured on a subboundary for different choices of u 0 . In the present paper, we improve the two above mentionned results by showing that the electric potential and the magnetic field can be stably and simultaneously recovered from a finite number of boundary observations of the solution. We stress out that the simultaneous identification of the magnetic field and the electric potential in the Schrödinger equation cannot be directly obtained from the results of [START_REF] Baudouin | Uniqueness and stability in an inverse problem for the Schrödinger equation[END_REF] and [START_REF] Cristofol | Stability estimate in an inverse problem for non autonomous magnetic Schrödinger equations[END_REF], as the electric (resp. magnetic) potential is a two (resp. first) order perturbation of the laplacian. As a matter of fact, the method of derivation of the stability estimate given in Theorem 1.2 is different for the one of [START_REF] Cristofol | Stability estimate in an inverse problem for non autonomous magnetic Schrödinger equations[END_REF][Theorem 1.1], as second order time-derivatives of the solution only are used.

1.3. Well posedness of the magnetic Schrödinger equation and main results. Before stating our main result, we need first to justify the existence of a unique solution of (1.1). To this end, we introduce the space

H 1 pΩq : H 1 0 pΩq equipped with the scalar product p¡∆ 1q 1{2 u, p¡∆ 1q 1{2 v ¡ L 2 pΩq , for any u, v H 1 pΩq,
and denote by H 2 : H 2 pΩq H 1 0 pΩq equipped with the scalar product p¡∆ 1qu, p¡∆ 1qv ¡ L 2 pΩq for any u, v H 2 pΩq.

Here and below ¤, ¤ ¡ L 2 pΩq denotes the usual scalar product in L 2 pΩq. Then, we have the following theorem:

Theorem 1.1. Let a A, q L V pΩ, Rq, χ C 3 p0, T ; Rq and β C 3 p0, T ; R). Then, for every u 0 satisfying ∆ k u 0 H 2 , k 0, 1, 2, and for any f W 2 p0, T ; H 2 pΩqq, there exists a unique solution u C 2 p0, T ; H 2 q C 3 p0, T ; H 0 q to the equation

(1.3) 4 p¡if t H a ptq βptqqpxqqupx, tq f, in Q Ω ¢ p0, T q, upx, T 2 q u 0 pxq in Ω.
Moreover, there exists a constant C ¡ 0 such that

(1.4) }f j t up., tq} H 1 pΩq ¤ C j ķ0 }∆ k u 0 } H 1 pΩq , j 0, 1, 2, t p0, T q.
Proof. Since H a ptq is a self adjoint operator in L 2 pΩq, associated with the sesquilinear form u Þ Ñ }pi∇ χptqaq} 2 L 2 pΩq n , u H 1 0 pΩq, then, it holds true (see. e. g. [START_REF] Cristofol | Stability estimate in an inverse problem for non autonomous magnetic Schrödinger equations[END_REF]) that the domain of H a ptq is DpH a ptqq H 1 0 pΩq H 2 pΩq. Further, as q L V pΩq for all t p0, T q, we deduce from the Kato-Rellich Theorem that DpH a ptq βptqqq DpH a ptqq. In view of [START_REF] Liu | Global uniqueness in determining electric potentials for a system of strongly coupled Schrödinger equations with magnetic potential terms[END_REF][Section ], there exists a family of unitary operotors pUpt, sqq 0¤s,t¤T in H 0 satisfying the following statements:

(1) U ps, sq Id, the identity mapping in H 0 , (2) U pt, sq DpH a psq βpsqqq DpH a ptq βptqqq, t, s r0, T s.

(3) For all φ DpH a psq βpsqqq, the mapping t Þ Ñ U pt, sqx, is continuously differentiable in r0, T s and satisfies ¡if t U pt, sq φ pH a ptq βptqqpxqqU pt, sqφ 0, therefore, arguing as in [START_REF] Cristofol | Stability estimate in an inverse problem for non autonomous magnetic Schrödinger equations[END_REF][Section 2], we check that up¤, tq U pt, T 2 qu 0 i

» t T 2 U pt, sqf psq ds, is solution to (1.
3) and satisfies the estimate 1.4. We complete the proof of this theorem..

In order to express the main result of this paper, let us introduce some notations. Let us denote by V, an arbitrary neighborhood of the boundary Γ. For M ¡ 0, and pa 0 , q 0 q A ¢L V pΩq, we define the admissible set of the unknown coefficients a and q:

S M pa 0 , q 0 q : tpa, qq A ¢ L V pΩq, such that a a 0 , and q q 0 in Vu.

Then, our main result can be stated as follows Theorem 1.2. Let M ¡ 0, let χ and β be as in Theorem 1.1 and satisfy 1.1. Let pa j , q j q, j 1, 2 be in S M pa 0 , q 0 q, where pa 0 , q 0 q are the same as above. Then, there exists n 1 initial conditions u 0,k , k 0, ..., n, such that we have

}a 1 ¡ a 2 } L 2 pΩq }q 1 ¡ q 2 } L 2 pΩq ¤ C ¡ n ķ0 }f ν f 2 t u 1,k ¡ f ν f 2 t u 2,k } 2 L 2 p0,T;L 2 pΓ qq © .
Here C ¡ 0 is a constant depending only on Ω, T, χ and β and u j,k , j 1, 2, is the solution of (1.1) where u 0,k is substituted for u 0 .

1.4. Comments.

Notice that as in [START_REF] Liu | Global uniqueness in determining electric potentials for a system of strongly coupled Schrödinger equations with magnetic potential terms[END_REF], we impose the initial condition up¤, T 2 q u 0 in (1.1) at t T 2 and not t 0. This allows us to "symmetrize" the solution u to (1.1) around t T 2 , and consequently to apply the Carleman estimate of Proposition 2.1 to u over p0, T q.

The assumption ∇ ¤ a 0 is purely technical and does not restrict the generality of Theorem 1.2. Indeed, it is well known that the magnetic potential is not meaningful in physics. The physical relevant quantity is the "two-form" da fa, which coïncides with the magnetic field curl a when n 3. Actually, given the "magnetic field" b, we can always choose a divergence free a such that we have b da. This amounts to substituting a ∇ψ for a, where ψ H 1 pΩq is solution to the system 4 ¡∆ψ ∇ ¤ a in Ω, ψ 0 in fΩ. As in [START_REF] Cristofol | Stability estimate in an inverse problem for non autonomous magnetic Schrödinger equations[END_REF], we enforce homogeneous Dirichlet-boundary conditions to (1.1). These homogeneous Dirichlet conditions impose that q be known in the vicinity V of the boundary fΩ. Nevertheless, this condition can be removed upon selecting suitable non-homogeneous Dirichlet boundary conditions on fΩ as in [START_REF] Baudouin | Uniqueness and stability in an inverse problem for the Schrödinger equation[END_REF]. Similarly, we can remove the assumption that a be known on V by selecting the initial conditions u 0,k , for k 1, ..., n, as in [START_REF] Cristofol | Stability estimate in an inverse problem for non autonomous magnetic Schrödinger equations[END_REF][Theorem 1.1]. Nevertheless, the set of initial conditions cannot be defined explicitly, so we rather stick with the formulation of Theorem1.2 given in this paper. Nevertheless, in order to avoid the inadequate expense of the site of this article, we shall not go further into details in this matter.

Evidently, it can be checked that if a 1 a 2 then, the electric potential can be Lipschitz stably retrieved from one boundary observation of the solution. This extends the result of Baudouin and Puel [START_REF] Baudouin | Uniqueness and stability in an inverse problem for the Schrödinger equation[END_REF], to the case of a magnetic Laplacian.

Similarly, if q 1 q 2 , we can determine the divergence free magnetic potential for n boundary observations of the solution, which generalises the result of [START_REF] Cristofol | Stability estimate in an inverse problem for non autonomous magnetic Schrödinger equations[END_REF].

It is worth mentionning that the stability estimate of Theorem 1.2 determines n 1 unknown functions pa 1 , ..., a n q and q from the knowledge of n 1 boundary observations over the time-space p0, T q. We stress out that the condition 1.1 imposed on the functions χ and β are essential in order to solve the inverse problem under study in this work. Indeed they allow us to recover the information on a and q from the knowledge of the "initial" condition of the second order derivative of the linearized system associated with 1.1 (see the second line in 3.11).

1.5. Outline of the paper. The proof of Theorem 1.2 is based on the celebrated Bugkheim. Klibanov method [START_REF] Bukhgeim | Uniqueness in the large of a class of multidimensional inverse problems[END_REF] which is by means of a Carleman estimate designed for Schrödinger equations. In Section 2, we state the global Carleman estimate we use in the derivation of the stability result. Section 3 contains the proof of Theorem 1.2.

GLOBAL CARLEMAN ESTIMATE

In this section, we recall the global Carleman inequality for vanishing solution on the boundary Σ, which can be found in [START_REF] Baudouin | Uniqueness and stability in an inverse problem for the Schrödinger equation[END_REF][Section3]( see also [START_REF] Cristofol | Stability estimate in an inverse problem for non autonomous magnetic Schrödinger equations[END_REF] ). This estimate is the main tool needed for the derivation of Theorem 1.2.

Given the Schrödinger operator (2.5)

L : if t ∆,
we define a function ψ C 4 pΩ, R q, satisfying the following conditions:

(i) |∇ψpxq| ¥ β ¡ 0, d x Ω.
(ii) ∇ψ ¤ ν 0 for all x ΓzΓ (iii) hΛ 1 ¡ 0, hε ¡ 0 such that for all ξ R n , and for all λ ¡ Λ 1 , we have λ|∇ψ ¤ ξ| 2 D 2 ψpξ, ξq ¥ ε|ξ| 2 , where D 2 ψ f 2 ψ fx i fx j ¨1¤i,j¤n and D 2 ψpξ, ξq denotes the C n -scalar product of D 2 r ψξ with ξ.

Notice that there are actual functions ψ verifying the above assumptions, such as x Þ Ñ |x ¡ x 0 | 2 , for an arbitrary x 0 R n zΩ and a subboundary Γ tx Γ, px ¡ x 0 q ¤ ν ¥ 0u. Furthermore, for λ ¡ 0 the following weight functions:

(2.6)

θpx, tq e λψpxq tpT ¡ tq , and ηpx, tq α ¡ e λψpxq tpT ¡ tq , where α ¡ }e λψ } L V pΩq . Finaly, we introduce the two operators P 1 and P 2 acting in C V 0 pQq I , as follows:

(2.7)

P 1 : if t ∆ s 2 |∇η| 2
, and P 2 : is f t η 2s∇η ¤ ∇ sp∆ηq, in such a way that P 1 P 2 e ¡sη Le sη . Then, we have the following result: Proposition 2.1. Assume that ψ and Γ satisfy the above conditions. Let η and θ be as in (2.6), and let P j , j 1, 2 be defined by (2.7). Then, there are two constants s 0 ¡ 0 and C ¡ 0, depending only on T , Ω and Γ , such that the estimate s}e ¡sη ∇u} 2

L 2 pQq s 3 }e ¡sη u} 2 L 2 pQq j 1,2 }P j e ¡sη u} 2 L 2 pQq ¤ C ¡ s}e ¡sη θ 1{2 pf ν ψq 1{2 f ν u} 2 L 2 pΣ q }e ¡sη Lu} 2 L 2 pQq © ,
holds for all s ¥ s 0 , and for any function u L 2 p0, T ; H 1 0 pΩqq such that Lu L 2 pQq and f ν u L 2 p0, T ; L 2 pΓ qq. Here Σ T stands for Γ ¢ p0, T q.

STABILITY ESTIMATE

In this section, we derive the stability estimate for a and q appearing in the magnetic Schrödinger equation (1.1) of Theorem 1.2. Here and henceforth the symbol " I " stands for the differentiation with respect to t. p¡if t H a j ptq βptqq j qu j 0, in Q, u j p¤, T 2 q u 0 pxq, in Ω,

u j 0, on Σ.
Then, u u 1 ¡ u 2 is a solution to the following boundary value problem (3.9)

6 8 7 p¡if t H a 1 βptqq 1 qu f, in Q, up¤, T 2 q 0 in Ω, u 0 on Σ,
where f χpa 1 ¡ a 2 q ¤ p¡2i∇ ¡ χpa 1 a 2 qqu 2 ¡ βpq 1 ¡ q 2 qu 2 . By differentiating (3.9), we get (3.10)

6 8 7 p¡if t H a 1 βptqq 1 qv g : f I ¡ H I a 1 u ¡ β I q 1 u, in Q, vp¤, T 2 q 0, in Ω, v 0, on Σ with v f t u. Thus, w f t v is a solution to (3.11) 6 8 7 p¡if t H a 1 ptq βptqq 1 qw h : f P ¡ 2pH I a 1 β I q 1 qv ¡ pH P a 1 β P q 1 qu, in Q, wp¤, T 2 q 2χ I p T 2 qpa 1 ¡ a 2 qpxq ¤ ∇u 0 ¡ iβ I p T 2 qpq 1 ¡ q 2 qpxqu 0 ,
in Ω, w 0, on Σ 3.2. Preliminary estimates. We start by stating a powerful tool introduced by A. L. Bughkeim and M. V. Klibanov in [START_REF] Bukhgeim | Uniqueness in the large of a class of multidimensional inverse problems[END_REF].

Lemma 3.1. For d t1, ..., nu, let η be given by (2.7). Then, there exists a positive constant κ ¡ 0, depending only on T , such that we have

» T 0 » Ω e ¡2sηpx,tq § § § » t T 2 ppξ, xq dξ § § § 2 dx dt ¤ κ s }e ¡s{η p} L 2 pQq d , for every p L 2 pQ T q d .
We turn now to establishing the coming statement with the aid of Proposition 2.1 and the above lemma. Lemma 3.2. There exists s 1 ¡ 0 such that for any s ¥ s 1 , we have the following estimate

s 3 }e ¡sη w} 2 L 2 pQq }P 1 e ¡sη w} 2 L 2 pQq ¤ C ¡ }e ¡sη pa 1 ¡ a 2 q} 2 L 2 pQq n }e ¡sη pq 1 ¡ q 2 q} 2 L 2 pQq s }e ¡sη θ 1{2 pf ν ψq 1{2 f ν w} 2 L 2 pΣ q © ,
where C is a positive constant independent of s.

Proof. By applying Proposition 2.1 to the solution w, we find a constant C ¡ 0 such that s 3 }e ¡sη w} 2 L 2 pQq }P 1 e ¡sη w} 2

L 2 pQq s}e ¡sη ∇w} 2 L 2 pQq ¤ C ¡ }e ¡sη L w} 2 L 2 pQq s}e ¡sη θ 1{2 pf ν ψq 1{2 f ν w} 2 L 2 pΣ q © , s ¡ s 0 , (3.13)
where Lwpx, tq p¡hpx, tq 2iχptqa 1 pxq ¤ ∇ χ 2 ptqa 2 1 pxq βptqq 1 pxqqwpx, tq. Here hpx, tq is given by the following identity hpx, tq ¡2pH I a 1

β I q 1 qv ¡ pH P a 1 ptq β P ptqq 1 qu f 1 pq 1 ¡ q 2 qpxq f 2 pa 1 ¡ a 2 qpxq, where f 1 px, tq β P u 2 2β I f t u 2 βf 2 t u 2 , and

f 2 px, tq χ P ¡ ¡ 2i∇ ¡ χpa 1 a 2 q ¡ 2χ I2 pa 1 a 2 q ¡ χχ P pa 1 a 2 © u 2 2χ I ¡ ¡ 2i∇ ¡ χpa 1 a 2 q ¡ 2χχ I pa 1 a 2 q © f t u 2 χ ¡ ¡ 2i∇ ¡ χpa 1 ¡ a 2 q © f 2 t u 2
In view of (1.4), we have f j C 0 pr0, T s; L V pΩqq for j 1, 2. Moreover, it is easy to see that H I a 1

β I q 1 and H P a 1 β P q 1 are bounded operators from L 2 p0, T ; H 1 pΩqq into L 2 pQq. Thus, there exists C ¡ 0, independent of s, such that we have

s 3 }e ¡sη w} 2 L 2 pQq }P 1 e ¡sη w} 2 L 2 pQq s}e ¡sη ∇w} 2 L 2 pQq ¤ C ¡ }e ¡sη pa 1 ¡ a 2 q} 2 L 2 pQq n }e ¡sη pq 1 ¡ q 2 q} 2 L 2 pQq s}e ¡sη θ 1{2 pf ν ψq 1{2 f ν w} 2 L 2 pΣ q ρu,v,w p}e ¡sη ρ} 2 L 2 pQ T q }e ¡sη ∇ρ} 2 L 2 pQq q © .
Therefore, since ∇ k up¤, tq

» t T 2 ∇ k vp¤, τ q dτ and ∇ k vp¤, tq » t T 2
∇ k wp¤, τ q dτ , for k 0, 1, we deduce from Lemma 3.1 that

s 3 }e ¡sη w} 2 L 2 pQq }P 1 e ¡sη w} 2 L 2 pQq s}e ¡sη ∇w} 2 L 2 pQq ¤ C ¡ }e ¡sη pa 1 ¡ a 2 q} 2 L 2 pQq n }e ¡sη pq 1 ¡ q 2 q} 2 L 2 pQq s}e ¡sη θ 1{2 pf ν ψq 1{2 f ν w} 2 L 2 pΣ q }e ¡sη w} 2 L 2 pQ T q }e ¡sη ∇w} L 2 pQq © ,
for any s ¥ s 0 . Thus, by taking s sufficiently large, we obtain

s 3 }e ¡sη w} 2 L 2 pQq }P 1 e ¡sη w} 2 L 2 pQq ¤ C ¡ }e ¡sη pa 1 ¡ a 2 q} 2 L 2 pQq }e ¡sη pq 1 ¡ q 2 q} 2 L 2 pQq s }e ¡sη θ 1{2 pf ν ψq 1{2 f ν w} 2 L 2 pΣ q © .
This completes the proof of the Lemma.

¥ s 2 }φp¤, T 2 q} 2 L 2 pΩq 4χ I p T 2 q 2 }e ¡sηp¤, T 2 q pa 1 ¡ a 2 q ¤ ∇u 0 } 2 L 2 pΩq n β I p T 2 q 2 }e ¡sηp¤, T 2 q pq 1 ¡ q 2 qu 0 } 2 L 2 pΩq (3.16)
¤ Cs ¡3{2 ¡ }e ¡sη pa 1 ¡ a 2 q} 2 L 2 pQq n }e ¡sη pq 1 ¡ q 2 q} 2 L 2 pQq s }e ¡sη θ 1{2 pf ν ψq 1{2 f ν w} 2

L 2 pΣ q © ,
Let us now choose the initial conditions u 0 as follows. Pick ω Ω such that ω ΩzV. Then, we choose u 0 C 6 0 pΩq such that u 0 pxq 1 for any x ω. Taking into account that q 1 ¡ q 2 and a 1 ¡ a 2 vanish in V and that ηpx, T 2 q ¤ ηpx, tq for all x Ω, we deduce from (3.16) that (3.18)

C 1 }e ¡sηp¤, T 2 q pq 1 ¡q 2 q} 2 L 2 pΩq ¤ Cs ¡3{2 ¡ }e ¡sηp¤, T 2 q pa 1 ¡a 2 q} 2 L 2 pΩq n }e ¡sηp¤, T 2 q pq 1 ¡q 2 q} 2 L 2 pΩq s }f ν w} 2

L 2 pΣ q © .
Here we used the fact that θe ¡2sη and f ν ψ are bounded on Σ . Next, we select n initial conditions u 0,k C 6 0 pΩq, for k 1, ..., n, such that u 0,k x k on ω. Then, we infer from (3.16) in a similar way that (3.19)

C 2 }e ¡sηp¤, T 2 q pa 1 ¡a 2 q k } 2 L 2 pΩq ¤ Cs ¡3{2 ¡ }e ¡sηp¤, T 2 q pa 1 ¡a 2 q} 2 L 2 pΩq n }e ¡sηp¤, T 2 q pq 1 ¡q 2 q} 2 L 2 pΩq s }f ν w} 2

L 2 pΣ q © ,
where pa 1 ¡ a 2 q k denotes the k th component of a 1 ¡ a 2 . Summing up (3.18) with (3.19) for k 1, .., n, we get that C 1 }e ¡sηp¤, T 2 q pq 1 ¡ q 2 q} 2 L 2 pΩq C 2 n j1 }e ¡sηp¤, T 2 q pa 1 ¡ a 2 q j } 2 L 2 pΩq ¤ C pn 1q s ¡3{2 ¡ }e ¡sηp¤,0q pa 1 ¡ a 2 q} 2 L 2 pΩq n }e ¡sηp¤, T 2 q pq 1 ¡ q 2 q} 2 L 2 pΩq s }f ν w} 2 L 2 pΣ q © , s ¥ s 2 Thus, there exists s 3 ¡ 0 such that for s ¥ s 3 , we have }e ¡sηp¤, T 2 q pq 1 ¡ q 2 q} 2 L 2 pΩq }e ¡sηp¤, T 2 q pa 1 ¡ a 2 q} 2 L 2 pΩq n ¤ Cs ¡1{2 }f ν f 2 t pu 1 ¡ u 2 q} 2 L 2 pΣ q . Finally, from the inequality e ¡sηp¤, T 2 q ¥ e ¡4s α¡1 T 2 ¡ 0, we get the desired result.

3. 1 .

 1 Linearization. Let u j , for j 1, 2, be solutions to (

  3.3. Proof of Theorem1.2. Let us now complete the stability estimate. Putting φpx, tq e ¡sηpx,tq wpx, tq and using the fact that φpx, 0q 0, we get

	Therefore, we get from the Cauchy-Schwarz inequality that
	}φp¤, T 2 q} 2 L 2 pΩq ¤ 2}P 1 φ} 2 L 2 pQq }φ} 2 L 2 pQq
	(3.15)	¤ s ¡3{2 ¡	s 3 }e ¡sη w} 2 L 2 pQq }P 1 e ¡sη w} 2 L 2 pQq	©	, s ¡ 0.
	Then, by Lemma 3.2, we obtain for all s
	}φp¤, T 2 q} 2 L 2 pΩq	0 » T 2	»	Ω	f t |φpx, tq| 2 dx dt 2	0 ¡ » T 2	»	Ω	© f t φpx, tqφpx, tq dx dt	.
	Hence, from the Green formula and (2.7) one can see that
	}φp¤, T 2 q} 2 L 2 pΩq 2		¡ » T 2 0 ¡ » T 2	» »
					0

Ω pif t ∆ s 2 |∇ηpx, tq| 2 qφpx, tqφpx, tq dx dt © 2 Ω P 1 φpx, tqφpx, tq dx dt © .