A new semismooth Newton method for solving finite-dimensional quasi-variational inequalities

In this paper, we consider the numerical method for solving finite-dimensional quasi-variational inequalities with both equality and inequality constraints. Firstly, we present a semismooth equation reformulation to the KKT system of a finite-dimensional quasi-variational inequality. Then we propose a semismooth Newton method to solve the equations and establish its global convergence. Finally, we report some numerical results to show the efficiency of the proposed method. Our method can obtain the solution to some problems that cannot be solved by the method proposed in (Facchinei et al. in Comput. Optim. Appl. 62:85–109, 2015). Besides, our method outperforms than the interior point method proposed in (Facchinei et al. in Math. Program. 144:369–412, 2014).


Introduction
We consider the finite-dimensional quasi-variational inequality QVI(K , F): Find a vector x * ∈ K(x * ) such that F x * T yx * ≥ 0, ∀y ∈ K x * , (1.1) where F : R n → R n is a point to point mapping and K : R n ⇒ R n is a point to set mapping with closed and convex images. Throughout the paper, we assume that F belongs to C 1 and for each x ∈ R n , the feasible set mapping K is given by K(x) y ∈ R n | g(y, x) ≤ 0, h(y, x) = 0 , (1.2) where g : R n × R n → R m 1 belongs to C 2 and g i (·, x) is convex on R n for each i = 1, 2, . . . , m 1 and for all x ∈ R n , h : R n × R n → R m 2 belongs to C 2 and h j (·, x) is affine on R n for each j = 1, 2, . . . , m 2 and for all x ∈ R n . When the set K(x) is independent of x, (1.1) reduces to the famous variational inequality (VI). For VI, we refer the reader to [13] and the references therein. QVI (1.1), which was first introduced by Bensoussan and Lions [2,3], has important applications in many fields such as generalized Nash games, mechanics, economics, statistics, transportation and biology; see for example [1,6,10,12] and the references therein. One interesting topic on QVI is to develop the efficient algorithms for the solution of QVI. Since QVI is nonsmooth and nonconvex, it is difficult to design effective methods for QVI, and by now, compared with VI, the numerical methods are still scarce. In this paper, we mainly focus on the numerical method based on the KKT conditions of QVI. This area attracts many people's attention and much progress has been made. In [12] an interior point approach was proposed to solve QVI and the convergence was established for several classes of interesting QVIs. Reference [8,9] developed a so called LP-Newton method and the method can be successfully applied to nonsmooth systems of equations with non-isolated solutions. Reference [21] developed an efficient regularized smoothing Newton-type algorithm for QVI. The proposed algorithm takes the advantage of newly introduced smoothing functions and a non-monotone line search strategy. [10] proposed a semismooth Newton method for QVI. They obtained global convergence and locally superlinear/quadratic convergence result for some important classes of quasi-variational inequality problems. The numerical results show that the method performs well.
There are many ways to compute a numerical solution of the nonlinear complementarity problems (NCP), such as linearized projected relaxation methods [13], the modulus-based matrix splitting method [24] and the penalty method [7,23,25]. In the past two decades, the nonsmooth-equation-based method has been thoroughly studied to solve NCP; see for example [5,[14][15][16][17][18][19] and the references therein. A common way to reformulate the complementarity system is to use the so called NCP-function. A function φ : For example, the famous Fischer-Burmeister (FB) function takes the form By the use of the NCP-function, nonlinear complementarity problem can be easily converted into a system of nonlinear equations. Most existing NCP-functions are generally nondifferentiable in the sense of Fréderivative but semismooth in the sense of Mifflin [20] and Qi and Sun [22]. In [17], the authors proposed a nonsmooth equation reformulation to the NCP. Their reformulation enjoys a nice property that it is continuous differentiable everywhere except at the solution. In this paper, we present a semismooth equation reformulation to the KKT system of a quasi-variational inequality and propose a semismooth Newton method to solve the equations. The paper is organized as follows. In the next section, we describe a semismooth equation reformulation to the KKT system of a quasi-variational inequality, present the semismooth Newton method and establish the global convergence for the method. In Sect. 3, we compare the proposed method with some other methods on problems list in [11].
In the following, we introduce some notations that will be used in this paper. For a continuously differentiable function F : R n → R n , we write JF(x) for the Jacobian of F at a point x ∈ R n , whereas ∇F(x) denotes the transposed Jacobian of F. Given a smooth mapping g : R n × R n → R m , (y, x) → g(y, x), ∇ y g(y, x) denotes the transpose of the partial Jacobian of g with respect to the y-variables. If F is locally Lipschitz continuous around x, then ∂F(x) denotes Clarke's generalized Jacobian of F at x. For a vector x ∈ R n and a subset I ⊂ {1, 2, . . . , n}, we write x I for the subvector consisting of the elements x i , i ∈ I. For a matrix A ∈ R n×n and two subsets I, J ⊂ {1, 2, . . . , n}, the symbol A IJ stands for the submatrix with entries a ij for i ∈ I, j ∈ J. The symbol diag(a 11 , a 22 , . . . , a nn ) stands for a diagonal matrix with diagonal elements a 11 , a 22 , . . . , a nn .

Semismooth equation reformulation and semismooth Newton method
Firstly, we give the following definition that will be used. Vd exists for any d ∈ R n , where ∂F(x) is the generalized Jacobian of F at x. F is strongly semismooth at x ∈ R n if for any d → 0 and any V ∈ ∂F(x + d), where F (x; d) denotes the directional derivative of F at x along the direction d.
A point x is called a KKT point of QVI (1.1) if there exist Lagrange multipliers λ ∈ R m 1 and ν ∈ R m 2 such that Similar to Theorem 1 of [12], we find that x * ∈ K(x * ) is a solution of (1.1) if there exist λ * ∈ R m 1 and ν * ∈ R m 2 such that (x * , λ * , ν * ) satisfies the KKT conditions (2.1). Moreover, if x * ∈ K(x * ) is a solution of (1.1) and some suitable constraint qualification holds at x * , then there exist λ * ∈ R m 1 and ν * ∈ R m 2 such that (x * , λ * , ν * ) satisfies the KKT conditions (2.1). Based on the above relationship, our aim is to develop a numerical method for solving the KKT conditions (2.1). For convenience, let and then (2.1) can be rewritten as where the w ∈ R m 1 are slack variables. It is not easy to solve (2.2) directly since the fourth formula is a complementarity system. We replace the complementarity system by an NCP-function [17], which is called the smoothed form of FB function: It is clear that, for each ε = 0, φ(u, v, ε) is continuously differentiable. We use it to construct an almost smooth equation reformulation to the fourth formula.
is a solution of the nonsmooth system of equations Associated with the system of H(x, λ, ν, w) = 0, we consider its natural merit function where we set z := (x, λ, ν, w).
By a direct calculation, we find that the gradient ∇θ (λ, w) of θ (·, ·) at (λ, w) can be expressed as follows: If θ (λ, w) = 0, we can get by a direct calculation Otherwise, we have λ i ≥ 0, w i ≥ 0 and λ i w i = 0 for any i, which means that, if λ 2 i + w 2 i = 0, then Therefore, the partial generalized derivatives S(λ, w) can be expressed in the form of and where a 2 i + b 2 i + c 2 i ≤ 1, u ∈ R m 1 satisfies u = 1, E is a matrix whose elements are one, V λ and V w are diagonal matrices whose diagonal elements belong to ∂ λ i φ FB i (λ i , w i ) and ∂ w i φ FB i (λ i , w i ), respectively. On the basis of the above calculations, we have the following proposition. (2.4). Then the following statements hold:

Proposition 2.2 Let the mapping H be defined by
(a) If F is continuously differentiable and g, h are twice continuously differentiable, then H is semismooth and where U λ , U w is defined by (2.7) and (2.8), respectively. (b) If, in addition, JF, ∇ 2 g i (i = 1, . . . , m 1 ) and ∇ 2 h j (j = 1, . . . , m 2 ) are locally Lipschitz, then H is strongly semismooth. (c) Let the merit function be defined by (2.5). If F is continuously differentiable and g, h are twice continuously differentiable, then is continuously differentiable, and its gradient is given by for an arbitrary element V ∈ ∂H(z).
That is, there are no equality constraints in QVI (1.1). Similarly, we can formulate the above problem in terms of the nonsmooth system of equations whereL(x, λ) := F(x) + ∇ y g(x, x)λ. Similar to the Proposition 2.2, if F is continuously differentiable and g is twice continuously differentiable, thenH is semismooth and where U λ and U w are the same as in Proposition 2.2. Now, we present the semismooth Newton method for (1.1).
Step 2. Choose an arbitrary element V k ∈ ∂H(z k ), and compute d k as a solution of the linear system of equations (2.11) If either this system is not solvable or the sufficient decrease condition is not satisfied, then take d k := -∇ (z k ).
Step 3. Compute a stepsize t k as the maximum of the numbers β l k , l k = 0, 1, 2, . . . , such that the following Armijo condition holds: Step 4. Set z k+1 := z k + t k d k , k ← k + 1, and go to Step 1.

End.
Below, we establish the following global convergence theorem for Algorithm 1. In the following, we suppose the direction is always given by (2.11). Suppose {z k } → z * and ∇ (z * ) = 0, by (2.11), we have Noting that V k cannot be 0, otherwise H(z k ) = 0 and z k would be a stationary point. Hence, we have  On the other hand, by (2.12), we have ∇ (z k ) Td ≤ -ρ d p < 0, which contradicts (2.17).
Remark 2.2 The method proposed in [10] only considers the case of inequality constraints, while our method can solve QVI with both equality and inequality constraints.
Besides, as we will see in the next section, our method can solve some problems in QVILIB [11], which cannot be solved by the method proposed by [10].

Numerical experiments
In this section, we report the results obtained by Algorithm 1 on problems list in QVILIB. All the computations in this paper were done using Matlab 2014a on a computer with 8.00 GB RAM and 2.5 GHz CPU. We solved all 55 test problems whose detailed description can be found in [11]. For each problem we list -the x-part of the starting point (the number reported is the value of all components of the x-part of the starting point); -the number of iterations; -the number of evaluations of ; -the value of Y (x, λ, ν) at the termination. In order to perform the linear algebra involved, we used Matlab's linear system solver mldivide. If any entry of the solution given by mldivide is a NaN or it is equal to ±∞ or the sufficient decrease condition is not satisfied, then an anti gradient direction is used. We take μ = 10 -5 , β = 0.5, ρ = 10 -10 , σ = 0.01 and p = 2.1. We choose λ 0 = 0, ν 0 = 0 and w 0 = 0 for all problems. For (2.6), we choose a i = b i = c i = 0 when (λ i , w i ) = (0, 0) and θ = 0. Our aim is mainly to verify the reliability of the method, and compare the iteration numbers   with the results presented in [10]. In order to perform a fair computation with the results in [10], we choose the same stopping criterion, i.e., let and choose the termination criterion to be Y (x k , λ k , ν k ) ≤ 10 -4 . The iteration is also stopped if the number of iterations exceeds 500 or the stepsize t k computed at Step 3 is less than 10 -6 . We denote Algorithm 2.2 proposed in [10] by SSN, and compare our method with SSN. The results are list in Table 1. From Table 1, for problems that can be solved by SSN, they can also be solved by our method with almost the same iteration numbers except problems Box2B, Box3A, KunR11, KunR12, KunR21 and KunR22. However, our method can solve the problems RHS1A1, RHS1B1, RHS2A1, RHS2B1 and Wal3, which cannot be solved by SSN.
We also compare our method with the interior point method (denoted by IP) proposed in [12] from the iteration number and CPU time. For IP, we use the same parameters presented in [12] and the results are list in Table 2. From Table 2, we can see that our method is much more effective than IP for most problems.
We also consider other problems in QVILIB which are not test in Tables 1 and 2, including the QVIs with equality constraints, that is, Problems LunSS1 to Scrim12 in Table 3. As we can see from the table, Algorithm 1 can solve over half of those problems effectively.
We tried to make some modifications to the algorithm for those problems that cannot be solved. Specifically, when calculating the Jacobian ofH(z), we use JF(x) to approximate  JL. For now, we cannot prove the convergence of the modified algorithm. However, it is interesting to find that the modified algorithm can find a solution for some problems, such as MoveSet3A1, MoveSet3A2, MoveSet3B1 and MoveSet3B2. The results are presented in Table 4. Figures 1-4 display the performance of our method on the problems BiLin1A, Movset 4A1, OutZ40 and Wal2. The vertical axis in those figures represents the value of Y and the horizontal axis represents the iteration number. As we can see from the figures, with the increase of the iteration numbers, the value of Y decrease.

Conclusion Remarks
In this paper, we have studied the numerical solution of QVI. We obtain the KKT system of a QVI and present a semismooth Newton method to solve the equations. We also establish its global convergence. Numerical results show that the performance of the proposed algorithm is promising.