On a Nonsmooth Gauss–Newton Algorithms for Solving Nonlinear Complementarity Problems

Abstract

In this paper, we propose a new version of the generalized damped Gauss–Newton method for solving nonlinear complementarity problems based on the transformation to the nonsmooth equation, which is equivalent to some unconstrained optimization problem. The B-differential plays the role of the derivative. We present two types of algorithms (usual and inexact), which have superlinear and global convergence for semismooth cases. These results can be applied to efficiently find all solutions of the nonlinear complementarity problems under some mild assumptions. The results of the numerical tests are attached as a complement of the theoretical considerations.

1. Introduction

Let $F:R^n\to R^n$ and let $F_i$, $i=1,\dots ,n$ denote the components of F. The nonlinear complementarity problem (NCP) is to find $\mathit{x}\in R^n$ such that

$$\mathit{x}\ge 0,F\left(\mathit{x}\right)\ge 0\mathrm{and}{\mathit{x}}^{T}F\left(\mathit{x}\right)=0.$$

(1)

The ith component of a vector $\mathit{x}$ is represented by $x_i$. Solving (1) is equivalent to solving a nonlinear equation $G\left(\mathit{x}\right)=0$, where the operator $G:R^n\to R^n$ is defined by

$$G\left(\mathit{x}\right)=\left[\begin{array}{c}\phi (x_1,F_1\left(\mathit{x}\right))\\ \dots \\ \phi (x_n,F_n\left(\mathit{x}\right))\end{array}\right]$$

with some special function $\phi$. Function $\phi$ may have one of the following forms:

$$\begin{array}{ccc}\hfill {\phi }_1(a,b)& =& min\{a,b\};\hfill \\ \hfill {\phi }_2(a,b)& =& \sqrt{a^2+b^2}-a-b;\hfill \\ \hfill {\phi }_3(a,b)& =& \theta \left(\right|a-b\left|\right)-\theta \left(a\right)-\theta \left(b\right),\hfill \end{array}$$

where $\theta :R\to R$ is any strictly increasing function with $\theta \left(0\right)=0$, see [1].

The (NCP) problem is one of the fundamental problems of mathematical programming, operations research, economic equilibrium models, and in engineering sciences. A lot of interesting and important applications can be found in the papers of Harker and Pang [2] and Ferris and Pang [3]. We can find the most essential applications in:

• engineering—optimal control problems, contact or structural mechanics problems, structural design problems, or traffic equilibrium problems,
• equilibrium modeling—general equilibrium (in production or consumption), invariant capital stock, or game-theoretic models.

We borrow a technique used in solving some smooth problems. If g is a merit function of G, i.e., $g\left(\mathit{x}\right)=\frac{1}{2}G{\left(\mathit{x}\right)}^{T}G\left(\mathit{x}\right)$, then any stationary point of $g\left(\mathit{x}\right)$ is a least-squares solution of the equation $G\left(\mathit{x}\right)=0$. Then, algorithms for minimization are equivalent to algorithms for solving equations. The usual Gauss–Newton method (known also as the differential corrections method), presented by Ortega and Rheinboldt [4] in the smooth case, has the form

$${\mathit{x}}^{(k+1)}={\mathit{x}}^{\left(k\right)}-{\left[G^{\prime }{\left({\mathit{x}}^{\left(k\right)}\right)}^T{G}^{\prime }\left({\mathit{x}}^{\left(k\right)}\right)\right]}^{-1}{G}^{\prime }{\left({\mathit{x}}^{\left(k\right)}\right)}^{T}G\left({\mathit{x}}^{\left(k\right)}\right).$$

(2)

Local convergence properties of the Gauss–Newton method was discussed by Chen and Li [5], but only for some smooth case. The Levenberg–Marquardt method is also considered, which is a modified Gauss–Newton method, in some papers, e.g., [6] or [7]. Moreover, some comparison of semismooth algorithms for solving (NCP) problems has been made in [8].

In practice, we may also consider the damped Gauss–Newton method

$${\mathit{x}}^{(k+1)}={\mathit{x}}^{\left(k\right)}-{\omega }_k{\left[G^{\prime }{\left({\mathit{x}}^{\left(k\right)}\right)}^T{G}^{\prime }\left({\mathit{x}}^{\left(k\right)}\right)+{\lambda }_{k}I\right]}^{-1}{G}^{\prime }{\left({\mathit{x}}^{\left(k\right)}\right)}^{T}G\left({\mathit{x}}^{\left(k\right)}\right)$$

(3)

with parameters ${\omega }_k$ and ${\lambda }_k$. Parameter ${\omega }_k$ may be chosen to ensure suitable decrease of g. If ${\lambda }_k$ is positive for all k, then the inverse matrix in (3) always exists because $G^{\prime }{\left({\mathit{x}}^{\left(k\right)}\right)}^{T}G\left({\mathit{x}}^{\left(k\right)}\right)$ is a symmetric and positive semidefinite matrix. The method (3) has the important advantage: the search direction always exists, even if $G^{\prime }\left(\mathit{x}\right)$ is singular. Naturally, in the case of nonsmooth equations, some additional assumptions are needed to allow the use of some line search strategies and to ensure the global convergence. Because, in some cases, a function G is nondifferentiable, so the equation $G\left(\mathit{x}\right)=0$ will be nonsmooth, whereby the method (3) may be useless. Some version of the Gauss–Newton method for solving complementarity problems was also introduced by Xiu and Zhang [9] for generalized problems, but only for linear ones. Thus, for solving nonsmooth and nonlinear problems, we propose two new versions of a damped Gauss–Newton algorithm based on B-differential. The usual generalized method is a relevant extension of the work by Subramanian and Xiu [10] for a nonsmooth case. In turn, an inexact version is related to the traditional approach, which was widely studied, e.g., in [11]. In recent years, various versions of the Gauss–Newton method were discussed, although most frequently for solving nonlinear least-squares problems, e.g., in [12,13].

The paper is organized as follows: in the next section, we review some notions needed, such as B-differential, BD-regularity, semismoothness, etc. (Section 2.1). Next, we propose a new optimization problem-based methods for the NCP, transforming the NCP into an unconstrained minimization problem by employing a function ${\phi }_3$ (Section 2.2). We state its global convergence and superlinear convergence rate under appropriate conditions. In Section 3, we present the results of numerical tests.

2. Materials and Methods

2.1. Preliminaries

If F is Lipschitz continuous, the Rademacher’s theorem [14] implies that F is almost everywhere differentiable. Let the set of points, where F is differentiable, be denoted by $D_F$. Then, the B-differential (the Bouligand differential) of F at $\mathit{x}$ (introduced in [15]) is

$${\partial }_{B}F\left(\mathit{x}\right)=\left\{\underset{{\mathit{x}}^{\left(n\right)}\to \mathit{x}}{lim}{F}^{\prime }\left({\mathit{x}}^{\left(n\right)}\right),{\mathit{x}}^{\left(n\right)}\in D_F\right\},$$

where $F^{\prime }\left(\mathit{x}\right)$ denotes the usual Jacobian of F at $\mathit{x}$. The generalized Jacobian of F at $\mathit{x}$ in the sense of Clarke [14] is

$$\partial F\left(\mathit{x}\right)=conv{\partial }_{B}F\left(\mathit{x}\right)$$

We say that F is BD-regular at $\mathit{x}$, if F is locally Lipschitz at $\mathit{x}$ and if all $\mathit{V}\in {\partial }_{B}F\left(\mathit{x}\right)$ are nonsingular (regularity on account of B-differential). Qi proved (Lemma 2.6, [15]) that, if F is BD-regular at $\mathit{x}$, then a neighborhood N of $\mathit{x}$ and a constant $\mathrm{C}>0$ exist such that, for any $\mathit{y}\in N$ and $\mathit{V}\in {\partial }_{B}F\left(\mathit{y}\right)$, $\mathit{V}$ is nonsingular and

$$∥{\mathit{V}}^{-1}∥\le \mathrm{C}.$$

Throughout this paper, $∥·∥$ denotes the 2-norm.

The notion of semismoothness was originally introduced for functionals by Mifflin [16]. The following definition is taken from Qi and Sun [17]. A function F is semismooth at a point $x$, if F is locally Lipschitzian at $\mathit{x}$ and

$$\underset{\mathit{V}\in \partial F(\mathit{x}+\mathrm{t}{\mathit{h}}^{\prime }),{\mathit{h}}^{\prime }\to \mathit{h},\mathrm{t}\downarrow 0}{lim}\mathit{V}{\mathit{h}}^{\prime }$$

exists for any $\mathit{h}\in R^n$. F is also said semismooth at $\mathit{x}$, if it is directionally differentiable at $\mathit{x}$ and

$$\mathit{V}\mathit{h}-F^{\prime }\left(\mathit{x},\mathit{h}\right)=o\left(∥\mathit{h}∥\right).$$

Scalar products and sums of semismooth functions are still semismooth functions. Piecewise smooth functions and maximum of a finite number of smooth functions are also semismooth. The semismoothness is the almost usually seen assumption on F in papers dealing with nonsmooth equations because it implies some important properties for convergence analysis of methods in nonsmooth optimization.

If for any $\mathit{V}\in \partial F(\mathit{x}+\mathit{h})$, as $\mathit{h}\to 0$

$$\mathit{V}\mathit{h}-F^{\prime }\left(\mathit{x},\mathit{h}\right)=O\left({∥\mathit{h}∥}^{1+p}\right),$$

where $0<\mathrm{p}\le 1$, then we say F is $\mathrm{p}$-order semismooth at $x$. Clearly, $\mathrm{p}$-order semismoothness implies semismoothness. If $\mathrm{p}=1$, then the function F is called strongly semismooth. Piecewise $C^2$ functions are examples of strongly semismooth functions.

Qi and Sun [17] remarked that, if F is semismooth at $\mathit{x}$, then, for any $\mathit{h}\to 0$

$$F(\mathit{x}+\mathit{h})-F\left(\mathit{x}\right)-F^{\prime }(\mathit{x};\mathit{h})=o\left(∥\mathit{h}∥\right),$$

and, if F is $\mathrm{p}$-order semismooth at $\mathit{x}$, then for any $\mathit{h}\to 0$

$$F(\mathit{x}+\mathit{h})-F\left(\mathit{x}\right)-F^{\prime }(\mathit{x};\mathit{h})=O\left({∥\mathit{h}∥}^{1+\mathrm{p}}\right).$$

Remark 1.

Strong semismoothness of the appropriate function usually implies quadratic convergence of method instead of the superlinear one for semismooth function.

In turn, Pang and Qi [18] proved that semismoothness of F at $\mathit{x}$ implies that

$$\underset{\mathit{V}\in \partial F(\mathit{x}+\mathit{h})}{sup}\left\{F(\mathit{x}+\mathit{h})-F\left(\mathit{x}\right)-\mathit{V}\mathit{h}\right\}=o\left(∥\mathit{h}∥\right).$$

Moreover, if F is $\mathrm{p}$-order semismooth at $\mathit{x}$, then

$$\underset{\mathit{V}\in \partial F(\mathit{x}+\mathit{h})}{sup}\left\{F(\mathit{x}+\mathit{h})-F\left(\mathit{x}\right)-\mathit{V}\mathit{h}\right\}=O\left({∥\mathit{h}∥}^{1+\mathrm{p}}\right).$$

2.2. The Algorithm and Its Convergence

Consider nonlinear equation $G\left(\mathit{x}\right)=0$ defined by ${\phi }_3$. The equivalence of solving this equation and problem (NCP) is described by the following theorem:

Theorem 1

(Mangasarian [1]). Let θ be any strictly increasing function from R into R, that is, $a>b\iff \theta \left(a\right)>\theta \left(b\right)$, and let $\theta \left(0\right)=0$. Then, $\mathit{x}$ solves the complementarity problem (1) if and only if

$$\theta \left(\left|F_i\left(\mathit{x}\right)-x_i\right|\right)-\theta \left(F_i\left(\mathit{x}\right)\right)-\theta \left(x_i\right)=0,i=1,2,\dots ,n.$$

(4)

For the convenience, denote

$$G_i\left(\mathit{x}\right):=\theta \left(\left|F_i\left(\mathit{x}\right)-x_i\right|\right)-\theta \left(F_i\left(\mathit{x}\right)\right)-\theta \left(x_i\right)$$

(5)

for $i=1,2,\dots ,n$.

We assume that the function $\theta$ in Theorem 1 has the form

$$\theta \left(\xi \right)=\xi \left|\xi \right|.$$

Let $G\left(\mathit{x}\right)$ be the associated function. We define function g in the following way:

$$g\left(\mathit{x}\right)=\frac{1}{2}{∥G\left(\mathit{x}\right)∥}^2,$$

which allows for solving system $G\left(\mathit{x}\right)=0$ based on solving the nonlinear least-square problem

$$\underset{\mathit{x}}{min}g\left(\mathit{x}\right).$$

(6)

Let us note that ${\mathit{x}}^{*}$ solves $G\left(\mathit{x}\right)=0$ if and only if it is a stationary point of g. Thus, from Theorem 1, ${\mathit{x}}^{*}$ solves (1).

Remark 2. 

On the other hand, the first-order optimality conditions for problem (6) are equivalent to the nonlinear system

$$\nabla g\left(\mathit{x}\right)=G^{\prime }{\left(\mathit{x}\right)}^{T}G\left(\mathit{x}\right)=0,$$

where $\nabla g$ is the gradient of g, provided G is differentiable and $G^{\prime }$ is the Jacobian matrix of G.

The continuous differentiability of the merit function g for some kind of nonsmooth functions was established by Ulbrich in the following lemma:

Lemma 1

(Ullbrich, [19]). Assume that the function $G:R^n\supset D\to R^n$ is semismooth, or, stronger, p-order semismooth, $0<p\le 1$, then the merit function $\frac{1}{2}{∥G\left(\mathit{x}\right)∥}^2$ is continuously differentiable on D with gradient $\nabla g\left(\mathit{x}\right)={\mathit{V}}^{T}G\left(\mathit{x}\right)$, where $\mathit{V}\in \partial G\left(\mathit{x}\right)$ is arbitrary.

Lemma 2.

For any $\mathit{x}\in R^n$, let ${\mathit{A}}_x={\mathit{V}}_x^T{\mathit{V}}_x$, where ${\mathit{V}}_x\in {\partial }_{B}G\left(\mathit{x}\right)$. Suppose that $\nabla g\left(\mathit{x}\right)\ne 0$. Then, given $\lambda >0$, the direction $\mathit{d}$ given by

$$({\mathit{A}}_x+\lambda \mathit{I})\mathit{d}=\nabla g\left(\mathit{x}\right)$$

is an ascent direction for g. In particular, there is a positive $\mathrm{w}$ such that $g(\mathit{x}-\mathrm{w}\mathit{d})<g\left(\mathit{x}\right)$.

Proof. 

There exist constants $\beta \ge 0$ and $\gamma >0$ such that

$$\beta {∥\mathit{h}∥}^2\le {\mathit{h}}^T{\mathit{A}}_x\mathit{h}\le \gamma {∥\mathit{h}∥}^2\mathrm{for}\mathrm{all}\mathit{h}\in R^n,$$

because ${\mathit{A}}_x$ defined as ${\mathit{V}}_x^T{\mathit{V}}_x$ is symmetric and positive semidefinite.

It follows that

$$(\beta +\lambda ){∥\mathit{h}∥}^2\le {\mathit{h}}^T({\mathit{A}}_x+\lambda \mathit{I})\mathit{h}\le (\gamma +\lambda ){∥\mathit{h}∥}^2\mathrm{for}\mathrm{all}\mathit{h}\in R^n.$$

Since $\nabla g\left(\mathit{x}\right)\ne 0$, $\mathit{d}\ne 0$. If we take $\mathit{h}=\mathit{d}$, we obtain

$${\mathit{d}}^T\nabla g\left(\mathit{x}\right)\ge (\beta +\lambda ){∥\mathit{d}∥}^2>0.$$

It follows that $\nabla g\left(\mathit{x}\right)\mathit{d}>0$ and that $\mathit{d}$ is a ascent direction for g (Section 8.2.1 in [4]).  □

Now, we present the generalized version of the damped Gauss–Newton method for solving the nonlinear complementarity problem.

Algorithm 1: The damped Gauss-Newton method for solving NCP

Let $\beta ,\delta \in (0,1)$ be given. Let ${\mathit{x}}^{\left(0\right)}$ be a starting point. Given ${\mathit{x}}^{\left(k\right)}$, the steps for obtaining ${\mathit{x}}^{(k+1)}$ are: Step 1: If $\nabla g\left({\mathit{x}}^{\left(k\right)}\right)=0$, then stop. Otherwise, choose any matrix ${\mathit{V}}_k\in {\partial }_{B}G\left({\mathit{x}}^{\left(k\right)}\right)$ and let ${\mathit{A}}_k={\mathit{V}}_k^T{\mathit{V}}_k$. Step 2: Let ${\lambda }_k=g\left({\mathit{x}}^{\left(k\right)}\right)$. Step 3: Find ${\mathit{d}}^{\left(k\right)}$ that is a solution of the linear system $$({\mathit{A}}_k+{\lambda }_k\mathit{I}){\mathit{d}}^{\left(k\right)}=\nabla g\left({\mathit{x}}^{\left(k\right)}\right).$$ Step 4: Compute the smallest nonnegative integer $m_k$ such that $$g({\mathit{x}}^{\left(k\right)}+{\beta }^m{\mathit{d}}^{\left(k\right)})-g\left({\mathit{x}}^{\left(k\right)}\right)\le -\delta {\beta }^m\nabla g{\left({\mathit{x}}^{\left(k\right)}\right)}^T{\mathit{d}}^{\left(k\right)}$$   and set $${\mathit{x}}^{(k+1)}={\mathit{x}}^{\left(k\right)}+{\beta }^{m_k}{\mathit{d}}^{\left(k\right)}.$$

Remark 3. 

(i) In Step 2, letting ${\lambda }_k=g\left({\mathit{x}}^{\left(k\right)}\right)$ is one of the simplest strategy because then $\left\{{\lambda }_k\right\}$ converges to 0.

(ii) The line search step (Step 4) in the algorithm follows the Armijo rule.

Theorem 2.

Let ${\mathit{x}}^{\left(0\right)}$ be a starting point and $\left\{x^{\left(k\right)}\right\}$ be a sequence generated by Algorithm 1. Assume that:

(a) ${sup}_k∥{\mathit{V}}_k∥<\infty$ for all ${\mathit{V}}_k\in {\partial }_{B}G\left({\mathit{x}}^{\left(k\right)}\right)$;

(b) $\nabla g\left(\mathit{x}\right)$ is Lipschitzian with a constant ${\mathrm{L}}_g>0$ on the level set $\mathrm{L}=\left\{\mathit{x}:g\left(\mathit{x}\right)\le g\left({\mathit{x}}^{\left(0\right)}\right)\right\}$.

Then, the generalized damped Gauss–Newton method described by Algorithm 1 is well defined and either $\left\{{\mathit{x}}^{\left(k\right)}\right\}$ terminates at a stationary point of g, or else every accumulation point of $\left\{{\mathit{x}}^{\left(k\right)}\right\}$, if it exists, is a stationary point of g.

Proof. 

The proof is almost the same as Theorem 2.1 in [10], providing appropriately modified assumptions.  □

For the nonsmooth case, the alternative condition may be considered instead of Lipschitz continuity of $\nabla g\left(\mathit{x}\right)$ (similar as in [10]). Thus, we have the following convergence theorem:

Theorem 3.

Let ${\mathit{x}}^{\left(0\right)}$ be a starting point and $\left\{x^{\left(k\right)}\right\}$ be a sequence generated by Algorithm 1. Assume that:

(a) the level set $\mathrm{L}=\left\{\mathit{x}:g\left(\mathit{x}\right)\le g\left({\mathit{x}}^{\left(0\right)}\right)\right\}$ is bounded;

(b) G is semismooth on $\mathrm{L}$.

Then, the generalized damped Gauss–Newton method described by Algorithm 1 is well defined and either $\left\{{\mathit{x}}^{\left(k\right)}\right\}$ terminates at a stationary point of g, or else every accumulation point of $\left\{{\mathit{x}}^{\left(k\right)}\right\}$, if it exists, is a stationary point of g.

Now, we take up the rate of convergence of the considered algorithm. The following theorem shows suitable conditions in various cases.

Theorem 4.

Suppose that ${\mathit{x}}^{*}$ is a solution of problem (1), G is semismooth, and G is BD-regular at ${\mathit{x}}^{*}$. Then, there exists a neighborhood $N_{*}$ of ${\mathit{x}}^{*}$ such that, if ${\mathit{x}}^{\left(0\right)}\in N_{*}$ and the sequence $\left\{{\mathit{x}}^{\left(k\right)}\right\}$ is generated by Algorithm 1, we have:

(i) ${\mathit{x}}^{\left(k\right)}\in N_{*}$ for all k and the sequence $\left\{{\mathit{x}}^{\left(k\right)}\right\}$ is linear convergent to $x^{*}$;

(ii) if $\delta <0.5$, then the convergence is at least superlinear;

(iii) If G is strongly semismooth, then the convergence is quadratic.

Proof. 

The proof of similar theorem given by Subramanian and Xiu [10] is based on three lemmas, which have the same assumptions as theorem. Now, we present these lemmas in versions adapted to our nonsmooth case.  □

Lemma 3.

Assume that ${\mathit{d}}_x$ is a solution of the equation

$$({\mathit{A}}_x+{\lambda }_x\mathit{I}){\mathit{d}}_{\mathit{x}}=\nabla g\left(\mathit{x}\right),$$

where

$${\lambda }_x=g\left(\mathit{x}\right)and{\mathit{A}}_x={\mathit{V}}_x^T{\mathit{V}}_x$$

for some matrix ${\mathit{V}}_x$ taken from ${\partial }_{B}G\left(\mathit{x}\right)$. Then, there is a neighborhood $D_1$ of ${\mathit{x}}^{*}$ such that, for all $\mathit{x}\in D_1,$

$$∥\mathit{x}-{\mathit{d}}_x-{\mathit{x}}^{*}∥=o\left(∥\mathit{x}-{\mathit{x}}^{*}∥\right).$$

Lemma 4.

There is a neighborhood $D_2$ of ${\mathit{x}}^{*}$ such that, for all $\mathit{x}\in D_2$,

(a) $g\left(\mathit{x}\right)=\frac{1}{2}{(\mathit{x}-{\mathit{x}}^{*})}^T{\mathit{A}}_{*}(\mathit{x}-{\mathit{x}}^{*})+o\left({∥\mathit{x}-{\mathit{x}}^{*}∥}^2\right)$,

(b) $g\left(\mathit{x}\right)=\frac{1}{2}{(\mathit{x}-{\mathit{x}}^{*})}^T{\mathit{A}}_x(\mathit{x}-{\mathit{x}}^{*})+o\left({∥\mathit{x}-{\mathit{x}}^{*}∥}^2\right)$.

Lemma 5.

Suppose that the conditions of Lemma 1 hold. Then, there is a neighborhood $D_3$ of ${\mathit{x}}^{*}$ such that, for all $\mathit{x}\in D_3$,

$$g(\mathit{x}-{\mathit{d}}_x)-g\left(\mathit{x}\right)+\frac{1}{2}\nabla g{\left(\mathit{x}\right)}^T{\mathit{d}}_x\le o\left({∥\mathit{x}-{\mathit{x}}^{*}∥}^2\right).$$

The proofs of Lemmas 5 and 6 are almost the same as in [10]; however, in the proof of Lemma 4, we have to take into account the semismoothness and to use Lemma 1 to obtain the desired result.

At the same time, in a similar way, we may show a suitable rate of convergence.

Now, we consider the inexact version of the considered method, which computes an approximate step, using the nonnegative sequence of forcing terms to control the level of accuracy.

For the above inexact version of the algorithm, we can state the analogous theorems which are equivalents of Theorems 2–4. Based on our previous results, the proof can be carried out almost in the same way as that of theorems for the ’exact’ version of the method. However, the condition (7), implied by inexactness given in Step 3 of Algorithm 2, has to be considered. Thus, we omit both theorems as proofs here.

Algorithm 2: The inexact version of the damped Gauss-Newton method for solving NCP

Let $\beta ,\delta ,\theta \in (0,1)$ and ${\eta }_k\in [0,1)$ for all k given. Let ${\mathit{x}}^{\left(0\right)}\in R^n$ be a starting point. Given ${\mathit{x}}^{\left(k\right)}$, the steps for obtaining ${\mathit{x}}^{(k+1)}$ are: Step 1: If $\nabla g\left({\mathit{x}}^{\left(k\right)}\right)=0$, then stop. Otherwise, choose any matrix ${\mathit{V}}_k\in {\partial }_{B}G\left({\mathit{x}}^{\left(k\right)}\right)$ and let ${\mathit{A}}_k={\mathit{V}}_k^T{\mathit{V}}_k$. Step 2: Let ${\lambda }_k=g\left({\mathit{x}}^{\left(k\right)}\right)$. Step 3: Find ${\mathit{d}}^{\left(k\right)}$ that is a solution of the linear system $$∥({\mathit{A}}_k+{\lambda }_k\mathit{I}){\mathit{d}}^{\left(k\right)}+\nabla g\left({\mathit{x}}^{\left(k\right)}\right)∥\le {\eta }_k∥\nabla g\left({\mathit{x}}^{\left(k\right)}\right)∥.$$ (7) Step 4: If $$∥\nabla g({\mathit{x}}^{\left(k\right)}+{\mathit{d}}^{\left(k\right)})∥\le \theta ∥\nabla g\left({\mathit{x}}^{\left(k\right)}\right)∥$$   then let $${\mathit{x}}^{(k+1)}={\mathit{x}}^{\left(k\right)}+{\mathit{d}}^{\left(k\right)}$$   and go to Step 1. Step 5: Compute the smallest nonnegative integer $m_k$ such that $$g({\mathit{x}}^{\left(k\right)}+{\beta }^m{\mathit{d}}^{\left(k\right)})-g\left({\mathit{x}}^{\left(k\right)}\right)\le -\delta {\beta }^m\nabla g{\left({\mathit{x}}^{\left(k\right)}\right)}^T{\mathit{d}}^{\left(k\right)}$$   and set $${\mathit{x}}^{(k+1)}={\mathit{x}}^{\left(k\right)}+{\beta }^{m_k}{\mathit{d}}^{\left(k\right)}$$   and go to Step 1.

3. Numerical Results

In this section, we present results of our numerical experiments, obtained by coding both algorithms in Code:Blocks. We use double precision on an Intel Core i7 3.2 GHz running under the Windows Server 2016 operating system. We applied the generalized damped Gauss–Newton method to solve three nonlinear complementarity problems. In the following examples: $N_1$ and $N_2$ denote the number of performed iterations to satisfy the stopping criterion $\left|{\mathit{x}}^{(k+1)}-{\mathit{x}}^{\left(k\right)}\right|<{10}^{-7}$, using Algorithms 1 and 2, respectively. The forcing terms in Algorithm 2 were chosen as follows: ${\eta }_k={\left(10k\right)}^{-1}$ for all k.

Example 1

(from Kojima and Shindo [20]). Let the function $F:R^4\to R^4$ have the form

$$\begin{array}{ccc}\hfill F^1\left(\mathit{x}\right)& =& 3x_1^2+2x_1{x}_2+2x_2^2+x_3+3x_4-6,\hfill \\ \hfill F^2\left(\mathit{x}\right)& =& 2x_1^2+x_1+x_2^2+10x_3+2x_4-2,\hfill \\ \hfill F^3\left(\mathit{x}\right)& =& 3x_1^2+x_1{x}_2+2x_2^2+2x_3+9x_4-9,\hfill \\ \hfill F^4\left(\mathit{x}\right)& =& x_1^2+3x_2^2+2x_3+3x_4-3.\hfill \end{array}$$

Problem (NCP) with the above function F has two solutions:

$${\mathit{x}}^{*}={(1,0,3,0)}^{T}and{\mathit{x}}^{**}={(\sqrt{6}/2,0,0,0.5)}^T$$

for which

$$F\left({\mathit{x}}^{*}\right)={(0,31,0,4)}^{T}andF\left({\mathit{x}}^{**}\right)={\left(0,2+\frac{\sqrt{6}}{2},0,0\right)}^T.$$

Thus, ${\mathit{x}}^{*}$ is a non-degenerate solution of (NCP) because

$$\mathrm{L}:=\left\{i:x_i^{*}=0,F^i\left({\mathit{x}}^{*}\right)=0\right\}=\varnothing ,$$

but ${\mathit{x}}^{**}$ is a degenerate solution.

Depending upon the starting point, we obtained the convergence iteration process to both solutions (see

Table 1. Results for Example 1.

x${}^{\left(0\right)}$	${\mathit{N}}_1$	${\mathit{N}}_2$	Solution
${(1,0,0,0)}^T$	9	11	${\mathit{x}}^{**}$
${(0,0,1,0)}^T$	failed	18	${\mathit{x}}^{**}$
${(0,0,0,1)}^T$	failed	failed	-
${(1,0,1,0)}^T$	7	10	${\mathit{x}}^{*}$
${(1,0,0,1)}^T$	7	9	${\mathit{x}}^{**}$
${(1,0,1,-5)}^T$	6	8	${\mathit{x}}^{**}$

or Figure 1).

Example 2.

Let function $F:R^2\to R^2$ be defined as follows:

$$F\left(\mathit{x}\right)=\left[\begin{array}{c}2x_1+x_2^2-6\\ -x_1^2+4x_1+\frac{1}{2}{x}_2-3\end{array}\right].$$

Then, problem (NCP) has two solutions:

-

non-degenerate

$${\mathit{x}}^{*}={(0,6)}^{T}forwhichF\left({\mathit{x}}^{*}\right)={(30,0)}^T$$

-

degenerate

$${\mathit{x}}^{**}={(3,0)}^{T}forwhichF\left({\mathit{x}}^{**}\right)={(0,0)}^T.$$

Similar to Example 1, we obtained the convergence iteration process for both solutions, depending on the starting point (see

Table 2. Results for Example 2.

x${}^{\left(0\right)}$	${\mathit{N}}_1$	${\mathit{N}}_2$	Solution
${(0,0)}^T$	5	7	${\mathit{x}}^{*}$
${(1,0)}^T$	1	2	${\mathit{x}}^{**}$
${(0,1)}^T$	4	7	${\mathit{x}}^{*}$
${(1,-1)}^T$	2	4	${\mathit{x}}^{**}$
${(-1,1)}^T$	4	7	${\mathit{x}}^{**}$
${(5,5)}^T$	3	5	${\mathit{x}}^{*}$
${(100,100)}^T$	3	6	${\mathit{x}}^{*}$

or Figure 2).

Example 3

(from Jiang and Qi [21]). Let function $F:R^n\to R^n$ has the form $F\left(\mathit{x}\right)=\mathit{M}\mathit{x}+\mathit{q}$, where

$$\mathit{M}=\left[\begin{array}{cccccc}4& -1& 0& \dots & 0& 0\\ -1& 4& -1& \dots & 0& 0\\ 0& -1& 4& \dots & 0& 0\\ \dots & \dots & \dots & \dots & \dots & \dots \\ 0& 0& 0& \dots & 4& -1\\ 0& 0& 0& \dots & -1& 4\end{array}\right],\mathit{q}={\left(-1,\dots ,-1\right)}^T.$$

Because F is strictly monotonic, the proper problem (NCP) has exactly one solution.

Calculations have been made for various n with one starting point ${\mathit{x}}^{\left(0\right)}={(0,\dots ,0)}^T$. For all tests, we obtain the same number of iterations $N_1=3$ and $N_2=4$.

4. Conclusions

We have given the nonsmooth version of the damped generalized Gauss–Newton method presented by Subramanian and Xu [10]. The generalized Newton algorithms related to the Gauss–Newton method are well-known important tools for solving nonsmooth equations, which arise from various nonlinear problems such as nonlinear complementarity or variational inequality. These algorithms are especially useful when the problem has many variables. We have proved that the sequences generated by the methods are superlinearly convergent under mild assumptions. Clearly, the semismoothness and BD-regularity are sufficient to obtain only a superlinear convergence of methods, while strong semismoothness even gives quadratic convergence. However, if function G is not semismooth or BD-regular or the gradient of g is not Lipschitzian, the Gauss–Newton methods may be useless.

The performance of both methods was evaluated in terms of the number of iterations required. The analysis of the numerical results seems to indicate that the methods are usually reliable for solving semismooth problems. The results show that the inexact approach can produce a noticeable slowdown by the number of iterations (compare $N_1$ and $N_2$ in Figure 1 and Figure 2). In turn, an important advantage is that the algorithms allow us to find various solutions to the problem (this can be observed in two examples: the first and second one). However, if there are many solutions of the problem, then the relationship between the starting point and the obtained solution may be unpredictable.

Clearly, traditional numerical algorithms aren’t the only method for solving the nonlinear complementarity problems, regardless of the degree of nonsmoothness. Except for the methods presented in the paper and mentioned in the Introduction, some computational intelligence algorithms can be used to solve (NCP) problems, i.a., monarch butterfly optimization (see [22,23]), the earthworm optimization algorithm (see [24]), the elephant herding optimization (see [25,26]), or the moth search algorithm (see [27,28]). All of these approaches are bio-inspired metaheuristic algorithms.