Inexact fixed-point iteration method for nonlinear complementarity problems

Abstract

Based on the modulus decomposition, the structured nonlinear complementarity problem is reformulated as an implicit fixed-point system of nonlinear equations. Distinguishing from some existing modulus-based matrix splitting methods, we present a flexible modulus-based inexact fixed-point iteration method for the resulting system, in which the subproblem can be solved approximately by a linear system-solver. The global convergence of the proposed method is established by assuming that the system matrix is positive definite. Some numerical comparisons are reported to illustrate the applicability of the proposed method, especially for large-scale problems.

Introduction

This work concentres on the following structured nonlinear complementarity problems (abbreviated as SNCP): Find $u\in {\mathbb{R}}^n$ such that

$$\begin{array}{r}u\ge 0,f(u):=Au+q+\varphi (u)\ge 0,\mathrm{and}{u}^{T}f(u)=0,\end{array}$$

(1)

where $A\in {\mathbb{R}}^{n\times n}$ is a large scale and sparse positive definite matrix, $q\in {\mathbb{R}}^n$ , and $\varphi (u):{\mathbb{R}}^n\to {\mathbb{R}}^n$ is a continuous nonlinear mapping satisfying

$$\{\begin{array}{c}\Vert \varphi (u)-\varphi (v)\Vert \le {\ell }_{\varphi }\Vert u-v\Vert ,\mathrm{\forall }u,v\in {\mathbb{R}}_{+}^n,\\ {\ell }_{\varphi }+{\sigma }_{A^N}^{max}<{\lambda }_1^S.\end{array}$$

(2)

Here ${\ell }_{\varphi }>0$ is a constant, ${\lambda }_1^S$ is the smallest eigenvalue of matrix $A^S:=\frac{1}{2}(A+A^T)$ , and ${\sigma }_{A^N}^{max}$ denotes the greatest singular value of matrix $A^N:=\frac{1}{2}(A-A^T)$ . We assume throughout this paper that the solution set of problem (1) is non-nonempty.

The SNCP with the form (1) has been broadly used in optimization theory and real-world applications consisting of nonlinear programming, financial analysis, economic equilibrium and engineering simulation.^1–3 For instance, the famous Karush–Kuhn–Tucker conditions of nonlinear programming, ⁴ the advection-diffusion control problems, ⁵ the American-style option pricing problems,^6,7 and the moving boundary problems. ⁸ Such problems can be modeled as finite-dimensional nonlinear complementarity problems NCPs of weak nonlinearity by resorting to some discretization schemes (e.g. 5-point difference method, time-state discretization, and grid method). To pursue the high accuracy of the model, scholars usually take a very fine discretization operation on the aforementioned problems, which leads to a very large-scale SNCP (at least millions).

For the standard NCP, the equation-based approaches are state-of-the-art methods. Founded on the NCP function, we can reformulate the NCP as a system of nonlinear equations. Thereby, some smoothing-based and semismooth Newton-type algorithms^9–17 are recommended for the resulting system. Also, due to the equivalence between the standard NCP and the variational inequality problem over ${\mathbb{R}}_{+}^n$ , some interesting methods for variational inequalities and maximal operators are also available, see Yuan, ¹⁸ Noor and Bnouhachem, ¹⁹ Bnouhachem and Yuan, ²⁰ He et al. ²¹ Particularly, for the solution of linear complementarity problem (LCP for short, i.e., $\varphi (u)=0$ in problem (1)), Murty ²² introduced a modulus iteration method, which aims to avoid projection in the implicit fixed-point iterations. This idea was further developed by Bai. ²³ The author designed a general framework of the modulus-based matrix splitting (MBMS for short) method for LCP, and proved its global convergence by assuming that the system matrix (i.e., $A$ in LCP) is a positive definite matrix or an $H_{+}$ -matrix. Liu et al. ²⁴ accelerated the MBMS method when the system matrix is an $H_{+}$ -matrix. Moreover, various generalizations, variations, and accelerations of the MBMS method can be found in Zheng and Yin, ²⁵ Zheng et al., ²⁶ Ke et al. ²⁷

Recently, the MBMS-type methods have been extended to a class of NCPs with weak nonlinearity. By assuming that $A$ is a positive definite matrix or an $H_{+}$ -matrix and $\varphi (\cdot )$ is strongly monotone, an MBMS iteration method for SNCP was developed by Xia and Li. ²⁸ Ma and Huang ²⁹ presented a modified MBMS method, and established its global convergence under the conditions that the system matrix $A$ has an $H$ -splitting and $\varphi (\cdot )$ is global Lipschitz over ${\mathbb{R}}^n$ . Xie et al. ³⁰ proposed a two-step MBMS method for solving SNCP, in which the nonlinear term $\varphi (\cdot )$ is assumed to be monotone and differentiable, and its Jacobi matrix is diagonal and bounded above. Among the mentioned MBMS-type methods, the lower triangular and the upper triangular parts of system matrix $A$ are frequently utilized to construct some matrix splitting strategies, which helps us to solve subproblems at each iteration by using the Gauss-Seidel-type iterative methods. However, the splitting of nearly half of the system matrix in each iteration may increase the total computational cost.

To utilize the structure of matrix $A$ for accelerating the MBMS-type methods, we propose a modulus-based inexact fixed-point (MBIFP) method for solving problem (1). The nonlinear term $\varphi (\cdot )$ considered in this paper is assumed to be ${\ell }_{\varphi }$ -Lipschitz over ${\mathbb{R}}_{+}^n$ . This assumption is weaker than the global Lipschitz continuity over ${\mathbb{R}}^n$ required in Ma and Huang ²⁹ as well as the smoothness and diagonal Jacobi matrix assumed in Xie et al. ³⁰ The main contributions of this paper are listed in the following:

We propose an inexact fixed-point method based on the modulus decomposition, and establish the global convergence under some mild assumptions.

The subproblem in the proposed method is a linear system of equations. If $A$ is positive definite, the subproblem can be solved inexactly by some efficient linear system solvers, such as the Cauchy–Barzilai-Borwein (CBB) method and the conjugate gradient method.

Without any extra acceleration techniques and precondition schemes, some numerical experiments on large-scale SNCP display that the proposed method outperforms some existing MBMS methods.

The rest of the paper is outlined as follows. The MBIFP iteration method for problem (1) is proposed in Section “The Modulus-Based Inexact Fixed-Point Method”. In Section “Convergence Analysis”, we establish the global convergence of the proposed method. To show the effectiveness of the proposed method, some numerical comparisons are reported in Section “Numerical Experiments”. Finally, in Section “Conclusions” we conclude this paper with some final remarks.

The modulus-based inexact fixed-point method

As usually, the Euclidean norm $\Vert x\Vert$ is used for vector $x\in {\mathbb{R}}^n$ and the spectral norm $\Vert A\Vert$ is used for matrix $A\in {\mathbb{R}}^{n\times n}$ . The identity matrix is denoted by $I$ . Let $[n]:=\{1,2,\dots ,n\}$ and $T$ be the transpose operator of vector or matrix.

Let $x\in {\mathbb{R}}^n$ and

$$\{\begin{array}{cc}u& =|x|+x,\\ f(u)& =\mathrm{\Omega }(|x|-x),\end{array}$$

(3)

where $\mathrm{\Omega }$ is a positive scalar matrix, i.e., $\mathrm{\Omega }=\omega I$ with $\omega >0$ , and the absolute value $|\cdot |$ is component-wise. Combining the nonlinear term $f(u)$ in problem (1) and the system (3), we obtain from some simple algebra manipulations that

$$(\mathrm{\Omega }+A)x=(\mathrm{\Omega }-A)|x|-q-\varphi (|x|+x).$$

(4)

We next show the equivalence between the system (4) and the problem (1).

Theorem 1

If $x\in {\mathbb{R}}^n$ is a solution of system (4), then $u:=|x|+x$ is a solution of problem (1). Conversely, if $u\in {\mathbb{R}}_{+}^n$ is a solution of problem (1), then $x:={\displaystyle \frac{1}{2}(u-\frac{1}{\omega }f(u))}$ is a solution of system (4).

Proof

Suppose that $x\in {\mathbb{R}}^n$ is a solution of system (4). Then,

$$\mathrm{\Omega }(|x|-x)=A(|x|+x)+q+\varphi (|x|+x)=f(|x|+x).$$

By $\mathrm{\Omega }=\omega I$ ( $\omega >0$ ) and the fact that $0\le (|x|-x)\perp (|x|+x)\ge 0$ , we obtain

$$0\le f(|x|+x)\perp (|x|+x)\ge 0,$$

which implies that $u=|x|+x$ solves problem (1).

Conversely, suppose that $u\in {\mathbb{R}}_{+}^n$ is a solution of problem (1). Let $x:={\displaystyle \frac{1}{2}(u-\frac{1}{\omega }f(u))}$ . We have

If $u_i=0$ , $f_i(u)>0$ for all $i\in [n]$ , one has

$$(|x|+x{)}_i=\frac{f_i(u)}{2\omega }+\frac{-f_i(u)}{2\omega }=0=u_i,$$

and

$$(\mathrm{\Omega }(|x|-x){)}_i=\omega (\frac{f_i(u)}{2\omega }-\frac{-f_i(u)}{2\omega })=f_i(u).$$

If $u_i>0$ , $f_i(u)=0$ for all $i\in [n]$ , one has

$$(|x|+x{)}_i=\frac{u_i}{2}+\frac{u_i}{2}=u_i,$$

and

$$(\mathrm{\Omega }(|x|-x){)}_i=\omega (\frac{u_i}{2}-\frac{u_i}{2})=0=f_i(u).$$

If $u_i=0$ , $f_i(u)=0$ for all $i\in [n]$ , one has

$$(|x|+x{)}_i=0=u_i\text{and}(\mathrm{\Omega }(|x|-x){)}_i=0=f_i(u).$$

In summary, we have that $|x|+x=u$ and $\mathrm{\Omega }(|x|-x)=f(u)$ , which deduces that $f(|x|+x)=\mathrm{\Omega }(|x|-x)$ . Hence, $x={\displaystyle \frac{1}{2}(u-\frac{1}{\omega }f(u))}$ solves equation (4).    □

Based on the Theorem 1, we are ready to propose the MBIFP method for solving the system (4), and equivalently solving the problem (1). Essentially, the method can be treated as a generalized framework of the MBMS method introduced by Bai. ²³ Algorithm 1 below outlines our method in detail.

Algorithm 1.

Modulus based inexact fixed point method, MBIFP

s0.	Give a tolerance $\epsilon >0$ , an initial point $x_0\in {\mathbb{R}}^n$ and a nonnegative sequence $\{{\epsilon }_k\}\to 0$ . Set $k:=0$ , $u_k=|x_k|+x_k$ , and $\text{res}(x_k):=\Vert min\{u_k,f(u_k)\}\Vert$ .
s1.	Adopt a linear system solver (denoted by $\mathcal{L}(B,b(x_k))$ ) for linear system $Bx=b(x_k)$ . If $\Vert r_k\Vert \le {\epsilon }_k$ where $r_k=B{x}_{k+1}-b(x_k)$ , then accept the new iterate $$x_{k+1}:=\mathcal{L}(B,b(x_k)).$$ (5)
s2.	Let $u_{k+1}=|x_{k+1}|+x_{k+1}$ . If $\text{res}(x_{k+1})\le \epsilon$ , then stop and accept $u_{k+1}$ as an approximate solution of problem (1); otherwise, let $k:=k+1$ , and go to s1.

Remark 1

The system $Bx=b(x)$ represents an equivalent abstract formulation of the system (4). Let ${\epsilon }_k=0$ and take a splitting scheme on matrix $\mathrm{\Omega }+A$ , we shall obtain a corresponding solution method for problem (1). Specially, let $D$ , $U$ and $L$ be the diagonal part, the strictly upper triangular part, and the strictly lower triangular part of matrix $\mathrm{\Omega }+A$ , respectively.

If $B=D$ and $b(x)=(\mathrm{\Omega }-A)|x|-q-\varphi (|x|+x)-(U+L)x$ , then Algorithm 1 becomes the modulus-based Jacobi method.

If $B=D+L$ and $b(x)=(\mathrm{\Omega }-A)|x|-q-\varphi (|x|+x)-Ux$ , then Algorithm 1 reduces to the modulus-based Gauss-Seidel method.

If $B=D+\beta L$ and $b(x)=\beta [(\mathrm{\Omega }-A)|x|-q-\varphi (|x|+x)]+[(1-\beta )D-\beta U]x$ , where $\beta >1$ is a relaxation parameter, then Algorithm 1 denotes the modulus-based successive over-relaxation method.

There are many efficient inexact methods for the linear system $Bx=b(x_k)$ , such as Krylov subspace methods, projection-based methods, preconditioning-based methods, inexact Newton-type methods, and multi-grid methods, see Saad ³¹ for more details. In what follows, we shall discuss the two structures (i.e., symmetric and asymmetric) of system matrix $A$ , which will be useful to the convergence analysis of the proposed MBIFP method in the next section.

Case 1: Matrix $A$ is symmetric. If $A$ is symmetric and $\mathrm{\Omega }+A$ is symmetric positive definite, let $B=\mathrm{\Omega }+A$ and $b(x)=(\mathrm{\Omega }-A)|x|-q-\varphi (|x|+x)$ . Then, the linear system $Bx=b(x_k)$ has a unique solution.

Case 2: Matrix $A$ is asymmetric. If matrix $A$ is asymmetric, we split $A$ into

$$A=A^S+A^N,$$

where $A^S=\frac{A+A^T}{2}$ and $A^N=\frac{A-A^T}{2}$ are the symmetric and skew-symmetric parts of matrix $A$ , respectively. Then, system (4) is identical to

$$(\mathrm{\Omega }+A^S)x=(\mathrm{\Omega }-A)|x|-q-\varphi (|x|+x)-A^{N}x.$$

(6)

If $\mathrm{\Omega }+A^S$ is symmetric positive definite, let $B=\mathrm{\Omega }+A^S$ and $b(x)=(\mathrm{\Omega }-A)|x|-q-\varphi (|x|+x)-A^{N}x$ , we get a symmetric linear system $Bx=b(x_k)$ which also has a unique solution.

Convergence analysis

In this section, we prove that the sequence $\{u_k\}$ generated by the MBIFP method (Algorithm 1) converges to a solution of problem (1).

Proposition 1

Let $a>0$ , and $b_1\le b_2\le \dots \le b_n$ be positive real numbers. Then, it holds that

$$\underset{i\in [n]}{max}\frac{|a-b_i|}{a+b_i}=\{\begin{array}{cc}\frac{b_n-a}{a+b_n},& \mathrm{if}a\le \sqrt{b_1{b}_n}\\ \frac{a-b_1}{a+b_1},& \mathrm{if}a\ge \sqrt{b_1{b}_n}\end{array}$$

(7)

Proof.

According to the position relationship between the real number $a$ and the interval $[b_1,b_n]$ , we proved the assertion by cases.

Case 1: If $a\le b_1$ , since $a\le b_1\le b_2\le \dots \le b_n$ are positive real numbers, $\{{\displaystyle \frac{b_i-a}{b_i+a}\}}$ is a monotonously non-descending sequence. Hence ${\displaystyle \underset{i\in [n]}{max}\frac{|a-b_i|}{a+b_i}={\displaystyle \underset{i\in [n]}{max}\frac{b_i-a}{b_i+a}=\frac{b_n-a}{b_n+a}}}$ .

Case 2: If $a\ge b_n$ , then $\{{\displaystyle \frac{a-b_i}{a+b_i}\}}$ is monotonously non-increasing. Hence ${\displaystyle \underset{i\in [n]}{max}\frac{|a-b_i|}{a+b_i}=\underset{i\in [n]}{max}\frac{a-b_i}{a+b_i}=\frac{a-b_1}{a+b_1}}$ .

Case 3: If $a\in [b_1,b_n]$ . Let $b_m$ ( $m\in [n]$ ) be the largest real number such that $a\ge b_m$ . We have

$$\begin{array}{rl}\underset{i\in [n]}{max}\frac{|a-b_i|}{a+b_i}& =max\{\underset{i\in [n],i\le m}{max}\frac{|a-b_i|}{a+b_i},\underset{i\in [n],i>m}{max}\frac{|a-b_i|}{a+b_i}\}\\ & =max\{\underset{i\in [n],i\le m}{max}\frac{a-b_i}{a+b_i},\underset{i\in [n],i>m}{max}\frac{b_i-a}{b_i+a}\}\\ & =max\{\frac{a-b_1}{a+b_1},\frac{b_n-a}{b_n+a}\}\\ & =\{\begin{array}{cc}{\displaystyle \frac{b_n-a}{a+b_n},}& \mathrm{for}a\le \sqrt{b_1{b}_n},\\ {\displaystyle \frac{a-b_1}{a+b_1},}& \mathrm{for}a\ge \sqrt{b_1{b}_n}.\end{array}\end{array}$$

Hence, the assertion (6) is true.   □

Proposition 2

Suppose that the system matrix $A$ is positive definite, $\mathrm{\Omega }=\omega I$ and $\omega >0$ . The eigenvalues of matrix $A^S$ ordered non-decreasingly are denoted by ${\lambda }_1^S\le {\lambda }_2^S\le \dots \le {\lambda }_n^S$ . Let

$${\sigma }_{\mathrm{\Omega }}=\Vert {\left(A_{+}^S\right)}^{-1}(\mathrm{\Omega }-A^S)\Vert +2{\ell }_{\varphi }\Vert {\left(A_{+}^S\right)}^{-1}\Vert +2\Vert {\left(A_{+}^S\right)}^{-1}{A}^N\Vert ,$$

(8)

where $A_{+}^S:=\mathrm{\Omega }+A^S$ . If the condition (2) holds and

$$\omega >{\omega }_0:=\frac{1}{2}\sqrt{({\lambda }_1^S-\delta {)}^2+4\delta {\lambda }_n^S}-\frac{1}{2}({\lambda }_1^S-\delta )$$

(9)

where $\delta ={\ell }_{\varphi }+{\sigma }_{A^N}^{max}$ . Then, there exists a positive constant ${\sigma }_M>0$ such that ${\sigma }_{\mathrm{\Omega }}\le {\sigma }_M<1$ .

Proof

By (7), Proposition 1 and Cauchy–Schwarz inequality, we have that

$$\begin{array}{cc}{\sigma }_{\mathrm{\Omega }}& =\Vert {\left(A_{+}^S\right)}^{-1}(\mathrm{\Omega }-A^S)\Vert +2{\ell }_{\varphi }\Vert {\left(A_{+}^S\right)}^{-1}\Vert +2\Vert {\left(A_{+}^S\right)}^{-1}{A}^N\Vert \\ & \le {\displaystyle \underset{i\in [n]}{max}\frac{|\omega -{\lambda }_i^S|}{\omega +{\lambda }_i^S}+2\frac{{\ell }_{\varphi }}{\omega +{\lambda }_1^S}+2\frac{{\sigma }_{A^N}^{max}}{\omega +{\lambda }_1^S}}\\ & =\{\begin{array}{cc}{\displaystyle \frac{{\lambda }_n^S-\omega }{\omega +{\lambda }_n^S}+\frac{2{\ell }_{\varphi }+2{\sigma }_{A^N}^{max}}{\omega +{\lambda }_1^S},}& \mathrm{if}\omega \le \sqrt{{\lambda }_1^S{\lambda }_n^S},\\ {\displaystyle \frac{\omega -{\lambda }_1^S}{\omega +{\lambda }_1^S}+\frac{2{\ell }_{\varphi }+2{\sigma }_{A^N}^{max}}{\omega +{\lambda }_1^S},}& \mathrm{if}\omega \ge \sqrt{{\lambda }_1^S{\lambda }_n^S}.\end{array}\end{array}$$

Using condition (2) we get $\delta <{\lambda }_1^S$ . To straightforwardly solve the following inequalities

$$(c1)\{\begin{array}{c}\frac{{\lambda }_n^S-\omega }{\omega +{\lambda }_n^S}+\frac{2\delta }{\omega +{\lambda }_1^S}<1\\ \omega \le \sqrt{{\lambda }_1^S{\lambda }_n^S}\end{array},(c2)\{\begin{array}{c}\frac{\omega -{\lambda }_1^S}{\omega +{\lambda }_1^S}+\frac{2\delta }{\omega +{\lambda }_1^S}<1\\ \omega \ge \sqrt{{\lambda }_1^S{\lambda }_n^S}\end{array},$$

we get the desired results associated to $\omega$ as follows:

• For (c1) we have

  $$\frac{1}{2}\sqrt{({\lambda }_1^S-\delta {)}^2+4\delta {\lambda }_n^S}-\frac{1}{2}({\lambda }_1^S-\delta )<\omega \le \sqrt{{\lambda }_1^S{\lambda }_n^S}.$$

  (10)

Let ${\sigma }_M=\frac{{\lambda }_n^S-\omega }{\omega +{\lambda }_n^S}+\frac{2\delta }{\omega +{\lambda }_1^S}$ . By (c1) we have ${\sigma }_{\mathrm{\Omega }}\le {\sigma }_M<1$ .

For (c2) we have

$$\omega \ge \sqrt{{\lambda }_1^S{\lambda }_n^S}.$$

(11)

Let ${\sigma }_M=\frac{\omega -{\lambda }_1^S}{\omega +{\lambda }_1^S}+\frac{2\delta }{\omega +{\lambda }_1^S}$ . By (c2) we also have ${\sigma }_{\mathrm{\Omega }}\le {\sigma }_M<1$ .   □

Proposition 3

If the assumptions of Proposition 2 hold, then

$$\Vert B^{-1}\Vert <\frac{1}{{\omega }_0}.$$

(12)

Proof

By the assumptions, we have $\omega >{\omega }_0$ . Note for the spectral norm of matrix, we have $\Vert X\Vert ={\lambda }_{max}(X)$ the largest eigenvalue for symmetric and positive definite $X$ .

(1) If $A$ is symmetric and positive definite, then $B=A+\mathrm{\Omega }$ . We have

$$\Vert B^{-1}\Vert =\frac{1}{\Vert B\Vert }=\frac{1}{\Vert A+\omega I\Vert }<\frac{1}{\omega }<\frac{1}{{\omega }_0}.$$

(13)

(2) If $A$ is asymmetric, then $B=A^S+\mathrm{\Omega }$ . Since $A$ is positive definite, $A^S$ is symmetric positive definite. We get

$$\Vert B^{-1}\Vert =\frac{1}{\Vert B\Vert }=\frac{1}{\Vert A^S+\omega I\Vert }<\frac{1}{\omega }<\frac{1}{{\omega }_0}.$$

(14)

Combining (13) and (14) yields the desired result.   □

For representation convenience, we set $H(x)=B^{-1}b(x)$ . The following lemma is crucial for the global convergence analysis, which analyzes the Lipschitz continuity of $H(x)$ .

Lemma 1

Suppose that the assumptions of Propositions 2 hold. If $\varphi (u)$ is ${\ell }_{\varphi }$ -Lipschitz continuous on ${\mathbb{R}}_{+}^n$ , and $\omega >{\omega }_0$ , then $H(x)$ is $\ell$ -Lipschitz on ${\mathbb{R}}^n$ , where $\ell \in (0,1)$ .

Proof

Let $g(x):=|x|+x$ , then $g(x)\in {\mathbb{R}}_{+}^n$ and $\varphi (|x|+x)=\varphi (g(x))$ . By the first claim in condition (2) we get

$$\Vert \varphi (g(x))-\varphi (g(y))\Vert \le {\ell }_{\varphi }\Vert g(x)-g(y)\Vert \le 2{\ell }_{\varphi }\Vert x-y\Vert ,\mathrm{\forall }x,y\in {\mathbb{R}}^n.$$

(15)

Next, we shall prove the desired Lipschitz continuity of $H(x)$ from the following two cases.

(1) If $A$ is asymmetric, then

$$\begin{array}{rl}H(x)& =B^{-1}b(x)\\ & ={(\mathrm{\Omega }+A^S)}^{-1}[(\mathrm{\Omega }-A)|x|-q-\varphi (g(x))-A^{N}x]\\ & ={(\mathrm{\Omega }+A^S)}^{-1}[(\mathrm{\Omega }-A^S)|x|-q-\varphi (g(x))-A^{N}g(x)].\end{array}$$

By Proposition 2 and $\omega >{\omega }_0$ , there exists a constant $0<{\sigma }_M<1$ such that

$$\begin{array}{r}\Vert {(\mathrm{\Omega }+A^S)}^{-1}(\mathrm{\Omega }-A^S)\Vert +2{\ell }_{\varphi }\Vert {(\mathrm{\Omega }+A^S)}^{-1}\Vert \end{array}+2\Vert {(\mathrm{\Omega }+A^S)}^{-1}{A}^N\Vert \le {\sigma }_M.$$

(16)

Hence, for all $x,y\in {\mathbb{R}}^n$ , by (15) and (16), and using Cauchy–Schwarz inequality we obtain that

$$\begin{array}{rl}& \Vert H(x)-H(y)\Vert \\ & =\Vert {(\mathrm{\Omega }+A^S)}^{-1}[(\mathrm{\Omega }-A^S)|x|-q-\varphi (g(x))-A^{N}g(x)]\\ & -{(\mathrm{\Omega }+A^S)}^{-1}[(\mathrm{\Omega }-A^S)|y|-q-\varphi (g(y))-A^{N}g(y)]\Vert \\ & \le \Vert {(\mathrm{\Omega }+A^S)}^{-1}(\mathrm{\Omega }-A^S)\Vert |x|-|y|\Vert \\ & +\Vert {(\mathrm{\Omega }+A^S)}^{-1}{A}^N\Vert \Vert g(x)-g(y)\Vert \\ & +\Vert {(\mathrm{\Omega }+A^S)}^{-1}\Vert \Vert \varphi (g(x))-\varphi (g(y))\Vert \\ & \le \Vert {(\mathrm{\Omega }+A^S)}^{-1}(\mathrm{\Omega }-A^S)\Vert \Vert x-y\Vert \\ & +2\Vert {(\mathrm{\Omega }+A^S)}^{-1}{A}^N\Vert \Vert x-y\Vert \\ & +2{\ell }_{\varphi }\Vert {(\mathrm{\Omega }+A^S)}^{-1}\Vert \Vert x-y\Vert \le {\sigma }_M\Vert x-y\Vert .\end{array}$$

Thereby, $H(x)$ is Lipschitz on ${\mathbb{R}}^n$ with constant $\ell ={\sigma }_M<1$ .

(2) If $A$ is symmetric, then

$$H(x)=B^{-1}b(x)=(\mathrm{\Omega }+A{)}^{-1}[(\mathrm{\Omega }-A)|x|-q-\varphi (g(x))].$$

By the symmetry of $A$ , one can easily deduces that $A=A^S$ and $A^N=0$ . Hence, ${\lambda }_1={\lambda }_1^S$ and ${\lambda }_n={\lambda }_n^S$ . By Proposition 2 and $\omega >{\omega }_0$ , there exists a constant $0<{\sigma }_M<1$ satisfying

$$\begin{array}{rl}& \Vert {(\mathrm{\Omega }+A)}^{-1}(\mathrm{\Omega }-A)\Vert +2{\ell }_{\varphi }\Vert {(\mathrm{\Omega }+A)}^{-1}\Vert \\ & =\Vert {\left(\mathrm{\Omega }+A^S\right)}^{-1}\left(\mathrm{\Omega }-A^S\right)\Vert +2{\ell }_{\varphi }\Vert {\left(\mathrm{\Omega }+A^S\right)}^{-1}\Vert \\ & +2\Vert {\left(\mathrm{\Omega }+A^S\right)}^{-1}{A}^N\Vert \le {\sigma }_M.\end{array}$$

(17)

Then, for all $x,y\in {\mathbb{R}}^n$ , by (15), (17) and Cauchy–Schwarz inequality, we have

$$\begin{array}{rl}& \Vert H(x)-H(y)\Vert =\Vert {(\mathrm{\Omega }+A)}^{-1}[(\mathrm{\Omega }-A)|x|-q-\varphi (g(x))]\\ & -{(\mathrm{\Omega }+A)}^{-1}[(\mathrm{\Omega }-A)|y|-q-\varphi (g(y))]\Vert \\ & \le \Vert {(\mathrm{\Omega }+A)}^{-1}(\mathrm{\Omega }-A)\Vert |x|-|y|\Vert \\ & +\Vert {(\mathrm{\Omega }+A)}^{-1}\Vert \Vert \varphi (g(x))-\varphi (g(y))\Vert \\ & \le \Vert {(\mathrm{\Omega }+A)}^{-1}(\mathrm{\Omega }-A)\Vert \Vert x-y\Vert \\ & +2{\ell }_{\varphi }\Vert {(\mathrm{\Omega }+A)}^{-1}\Vert \Vert x-y\Vert \\ & \le {\sigma }_M\Vert x-y\Vert .\end{array}$$

Hence, $H(x)$ is also Lipschitz on ${\mathbb{R}}^n$ with constant $\ell ={\sigma }_M<1$ .   □

Lemma 2

Suppose that the assumptions of Propositions 2 hold, $\varphi (u)$ is ${\ell }_{\varphi }$ -Lipschitz on ${\mathbb{R}}_{+}^n$ , and sequence $\{x_k\}$ is generated by Algorithm 1 . Let $x^{*}=H(x^{*})$ and $\omega >{\omega }_0$ . Then we have

$$\Vert x_{k+1}-x^{*}\Vert \le {\ell }^{k+1}\Vert x_0-x^{*}\Vert +\frac{1}{{\omega }_0}\sum _{i=0}^k{\ell }^{k-i}{\epsilon }_i,\mathrm{\forall }k,$$

(18)

where $\ell \in (0,1)$ is the Lipschitz constant of $H(x)$ .

Proof.

By the step s1 in Algorithm 1, we have

$$x_{k+1}=B^{-1}b(x_k)+B^{-1}{r}_k=H(x_k)+B^{-1}{r}_k,\mathrm{\forall }k$$

where $\Vert r_k\Vert \le {\epsilon }_k$ . By Proposition 3, we obtain $\Vert x_{k+1}-H(x_k)\Vert <{\displaystyle \frac{{\epsilon }_k}{{\omega }_0}}$ . In what follows, we will prove (18) by induction. Let ${\eta }_k={\ell }^{k+1}\Vert x_0-x^{*}\Vert +{\displaystyle \frac{1}{{\omega }_0}\sum _{i=0}^k{\ell }^{k-i}{\epsilon }_i}$ .

For $k=0$ , by Lemma 1 and Proposition 3, we have

$$\begin{array}{rl}\Vert x_1-x^{*}\Vert & =\Vert x_1-H(x_0)+H(x_0)-H(x^{*})\Vert \\ & \le \Vert x_1-H(x_0)\Vert +\Vert H(x_0)-H(x^{*})\Vert \\ & \le {\displaystyle \frac{{\epsilon }_0}{{\omega }_0}+\ell \Vert x_0-x^{*}\Vert ={\eta }_0.}\end{array}$$

Assume that the formula (18) holds for $k\ge 1$ , it is $\Vert x_k-x^{*}\Vert \le {\eta }_{k-1}$ . Then, by Lemma 1 we have

$$\begin{array}{rl}\Vert x_{k+1}-x^{*}\Vert & =\Vert x_{k+1}-H(x_k)+H(x_k)-H(x^{*})\Vert \\ & \le \Vert x_{k+1}-H(x_k)\Vert +\Vert H(x_k)-H(x^{*})\Vert \\ & \le {\displaystyle \frac{{\epsilon }_k}{{\omega }_0}+\ell \Vert x_k-x^{*}\Vert \le {\displaystyle \frac{{\epsilon }_k}{{\omega }_0}+\ell {\eta }_{k-1}}}\\ & ={\displaystyle \frac{{\epsilon }_k}{{\omega }_0}+{\ell }^{k+1}\Vert x_0-x^{*}\Vert +\frac{1}{{\omega }_0}\sum _{i=0}^{k-1}{\ell }^{k-i}{\epsilon }_i}\\ & ={\ell }^{k+1}\Vert x_0-x^{*}\Vert +\frac{1}{{\omega }_0}\sum _{i=0}^k{\ell }^{k-i}{\epsilon }_i\\ & ={\eta }_{k+1},\end{array}$$

which completes the proof.   □

Lemma 3

If the assumptions of Lemma 2 hold, then

$$\begin{array}{r}\Vert x_{k+1}-x^{*}\Vert \le \frac{{\ell }^{k+1}}{1-\ell }(\Vert x_1-x_0\Vert +\frac{{\epsilon }_0}{{\omega }_0})+\frac{1}{{\omega }_0}\sum _{i=0}^k{\ell }^{k-i}{\epsilon }_i,\mathrm{\forall }k.\end{array}$$

(19)

Proof: By (18), it follows from Lemma 1 and Cauchy–Schwarz inequality that

$$\begin{array}{rl}\Vert x_0-x^{*}\Vert & =\Vert x_0-x_1+x_1-H(x_0)+H(x_0)-H(x^{*})\Vert \\ & \le \Vert x_1-x_0\Vert +\Vert x_1-H(x_0)\Vert +\Vert H(x_0)-H(x^{*})\Vert \\ & \le \Vert x_1-x_0\Vert +\frac{{\epsilon }_0}{{\omega }_0}+\ell \Vert x_0-x^{*}\Vert ,\end{array}$$

which implies that

$$\Vert x_0-x^{*}\Vert \le \frac{1}{1-\ell }(\Vert x_1-x_0\Vert +\frac{{\epsilon }_0}{{\omega }_0}).$$

(20)

Substituting (20) into the claim (18) yields the desired result (19).  □

The following theorem analyses the global convergence of Algorithm 1.

Theorem 2

Suppose that the assumptions of Lemma 2 hold, and $\{{\epsilon }_k\}\to 0$ . Then, the sequence $\{u_k=|x_k|+x_k\}$ generated by Algorithm 1 converges to a solution of problem (1).

Proof By Lemma 3, we have

$$\Vert x_{k+1}-x^{*}\Vert \le \frac{{\ell }^{k+1}}{1-\ell }(\Vert x_1-x_0\Vert +\frac{{\epsilon }_0}{{\omega }_0})+\frac{1}{{\omega }_0}\sum _{i=0}^k{\ell }^{k-i}{\epsilon }_i,\mathrm{\forall }k.$$

(21)

Where $0<\ell <1$ is the Lipschitz constant of $H(x)$ , and $x^{*}$ is a fixed point of $H(x)$ , i.e., $x^{*}=H(x^{*})$ , which is also a solution of problem (1).

Next, we only need to prove that the right-hand-side of inequality (21) converges to zero as $k\to \mathrm{\infty }$ . Since ${\epsilon }_k\to 0$ , for $\mathrm{\forall }ϵ>0$ sufficiently small, there exists a positive integer $K_1$ such that

$${\epsilon }_k<\frac{(1-\ell ){\omega }_0}{4}ϵ,\mathrm{\forall }k\ge K_1.$$

(22)

Let ${\epsilon }_{max}=max\{{\epsilon }_0,{\epsilon }_1,\dots ,{\epsilon }_{K_1-1}\}$ . We have from (22) that

$$\begin{array}{rl}\sum _{i=0}^k{\ell }^{k-i}{\epsilon }_i& ={\ell }^{k-(K_1-1)}\sum _{i=0}^{K_1-1}{\ell }^{K_1-1-i}{\epsilon }_i+\sum _{i=K_1}^k{\ell }^{k-i}{\epsilon }_i\\ & <{\ell }^{k-(K_1-1)}{\epsilon }_{max}\sum _{i=0}^{K_1-1}{\ell }^{K_1-1-i}+\frac{(1-\ell ){\omega }_0}{4}ϵ\sum _{i=K_1}^k{\ell }^{k-i}\\ & <{\ell }^{k-(K_1-1)}{\epsilon }_{max}\frac{1-{\ell }^{K_1-1}}{1-\ell }+\frac{(1-\ell ){\omega }_0}{4}ϵ\sum _{j=0}^{\mathrm{\infty }}{\ell }^j\\ & <\frac{{\ell }^{k-(K_1-1)}{\epsilon }_{max}}{1-\ell }+\frac{(1-\ell ){\omega }_0ϵ}{4}\frac{1}{1-\ell }\\ & =\frac{{\ell }^{k-(K_1-1)}{\epsilon }_{max}}{1-\ell }+\frac{{\omega }_0ϵ}{4}.\end{array}$$

Let ${\displaystyle \frac{{\ell }^{k-(K_1-1)}{\epsilon }_{max}}{1-\ell }<\frac{{\omega }_0ϵ}{4}}$ , we get $k>K_1-1+(\ln (1-\ell ){\omega }_0ϵ-\ln (4{\epsilon }_{max}))/\ln l$ . Let $K_2=\lceil (\ln (1-\ell ){\omega }_0ϵ-\ln (4{\epsilon }_{max}))/\ln l\rceil$ , then for all $k>K_1-1+K_2$ , we have ${\displaystyle \sum _{i=0}^k{\ell }^{k-i}{\epsilon }_i<\frac{{\omega }_0ϵ}{2}.}$ Since $0<\ell <1$ , there is positive integer $K_3$ such that for all $k>K_3$ ,

$$\frac{{\ell }^{k+1}}{1-\ell }(\Vert x_1-x_0\Vert +{\epsilon }_0)<\frac{ϵ}{2}.$$

Hence, for all $k\ge max\{K_1-1+K_2,K_3\}$ , we have that $\Vert x_{k+1}-x^{*}\Vert <ϵ.$ Let $ϵ\to 0$ , we have $\Vert x_{k+1}-x^{*}\Vert \to 0$ as $k\to \mathrm{\infty }$ and the assertion of Theorem 2 is proven.  □

Numerical experiments

In this section, some numerical experiments are presented for testing the effectiveness and performance of the proposed method. Specially, the CBB method proposed in Raydan and Svaiter ³² is adopted as a solver $\mathcal{L}(B,b(x_k))$ used in Algorithm 1. The CBB method is a combination of the Cauchy method ³³ and Barzilai–Borwein step-size. ³⁴ Algorithm 2 below summarizes the CBB solver $\mathcal{L}(B,b(x_k))$ in detail.

Algorithm 2.

Linear system solver $\mathcal{L}(B,b(x_k))$ : Cauchy-Barzilai-Borwein method, CBB

s0. Let $r(y)=By-b(x_k)$ . Set $y_0:=x_k$ , $m:=0$ .

s1.Let $r_m:=r(y_m)$ . If $\Vert r_m\Vert \le {\epsilon }_k$ , then stop and accept $x_{k+1}=y_m$ .

s2. Compute the BB stepsize by $$t_m=(r_m^T{r}_m)/(r_m^{T}B{r}_m).$$

s3. Update the new iterate by: $$y_{m+1}=z_m-t_{m}r(z_m),$$ where $z_m=y_m-t_m{r}_m$ .

s4. Let $m:=m+1$ , go to s1.

We denoted by MBIFP-CBB the Algorithm 1 cooperated with the CBB method (Algorithm 2). The MBIFP-CBB method is compared to the modulus-based successive over-relaxation (MSOR for short) method , two-step MSOR (TMSOR for short) method and two-step modulus-based accelerated over-relaxation (TMAOR for short) method. For details of the compared methods, see Xie et al. ³⁰ All the mentioned methods were coded by MATLAB R2016b and run on a PC with 1.60 GHz CPU and 8 GB RAM.

In the implementation of all methods, the stopping criterion is set to

$$\mathbf{r}\mathbf{e}\mathbf{s}:=\Vert min\{u_k,f(u_k)\}\Vert \le {10}^{-4},$$

and the initial point is given by $x_0=r\mathrm{and}(\mathrm{n},1)$ . The parameters are detailed in the following:

1. For the MSOR and TMSOR methods, three choices on parameter $\alpha$ are introduced in Xie et al., ³⁰ i.e., $\alpha =0.4$ , $\alpha =0.6$ , $\alpha =0.8$ denoted by MSOR ${}_a$ , MSOR ${}_b$ , MSOR ${}_c$ , and TMSOR ${}_a$ , TMSOR ${}_b$ , TMSOR ${}_c$ , respectively;

2. For the TMAOR method, two choices of $(\alpha ,\beta )$ are introduced in Xie et al., ³⁰ i.e., $\alpha =0.6$ , $\beta =0.2$ and $\alpha =0.2$ , $\beta =0.6$ , denoted by TMAOR ${}_a$ and TMAOR ${}_b$ , respectively;

3. For the proposed MBIFP-CBB method, we set ${\epsilon }_k=\frac{1}{k}$ . Figure 1 shows the values of ${\omega }_0=\frac{1}{2}\sqrt{({\lambda }_1^S-\delta {)}^2+4\delta {\lambda }_n^S}-\frac{1}{2}({\lambda }_1^S-\delta )$ associated with Examples 1 and 2, in which the scales vary from ${10}^2$ to ${100}^2$ . Based on the suggestive results in Figure 1, we set $\omega =3$ and $\omega =6$ for Example 1 and Example 2, respectively.

The notations used in Tables 1 and 2 are interpreted as follows: ind: performance indicators at termination of algorithm; n: dimension of test example; iter: total number of outer iterations; avg-init: average number of inner iterations; cput: CPU time in seconds; res: residual error defined by stopping criteria;

Figure 1.

The values of ${\omega }_0$ associated to Examples 1 and 2.

Table 1.

The numerical results on Example 1.

$ind$	$\mathit{n}$	MBIFP-CBB	MSOR ${}_a$	MSOR ${}_b$	MSOR ${}_c$	TMSOR ${}_a$	TMSOR ${}_b$	TMSOR ${}_c$	TMAOR ${}_a$	TMAOR ${}_b$
iter	${32}^2$	9	68	45	46	34	23	21	22	68
	${128}^2$	10	77	50	47	39	25	22	26	77
	${512}^2$	10	87	56	50	44	28	23	29	85
	${2048}^2$	11	96	62	53	48	31	25	32	94
avg-init	${32}^2$	1.1	1	1	1	2	2	2	2	2
	${128}^2$	1.4	1	1	1	2	2	2	2	2
	${512}^2$	1.4	1	1	1	2	2	2	2	2
	${2048}^2$	1.5	1	1	1	2	2	2	2	2
res	${32}^2$	5.7 $\times {10}^{-5}$	9.4 $\times {10}^{-5}$	8.1 $\times {10}^{-5}$	8.5 $\times {10}^{-5}$	9.6 $\times {10}^{-5}$	6.7 $\times {10}^{-5}$	9.4 $\times {10}^{-5}$	9.3 $\times {10}^{-5}$	9.5 $\times {10}^{-5}$
	${128}^2$	2.4 $\times {10}^{-5}$	9.9 $\times {10}^{-5}$	9.6 $\times {10}^{-5}$	9.8 $\times {10}^{-5}$	8.6 $\times {10}^{-5}$	9.6 $\times {10}^{-5}$	8.7 $\times {10}^{-5}$	6.8 $\times {10}^{-5}$	8.9 $\times {10}^{-5}$
	${512}^2$	9.6 $\times {10}^{-5}$	8.9 $\times {10}^{-5}$	9.1 $\times {10}^{-5}$	8.6 $\times {10}^{-5}$	7.7 $\times {10}^{-5}$	9.1 $\times {10}^{-5}$	9.5 $\times {10}^{-5}$	7.6 $\times {10}^{-5}$	9.9 $\times {10}^{-5}$
	${2048}^2$	8.0 $\times {10}^{-5}$	9.3 $\times {10}^{-5}$	8.7 $\times {10}^{-5}$	8.3 $\times {10}^{-5}$	9.3 $\times {10}^{-5}$	8.7 $\times {10}^{-5}$	7.1 $\times {10}^{-5}$	8.6 $\times {10}^{-5}$	9.6 $\times {10}^{-5}$
cput	${32}^2$	0.0007	0.0031	0.0025	0.0022	0.0025	0.0021	0.0017	0.0022	0.0077
	${128}^2$	0.0131	0.0726	0.0580	0.0610	0.1079	0.0603	0.0575	0.0693	0.2088
	${512}^2$	0.4852	2.8447	1.8949	1.5946	2.6805	1.8892	1.3160	1.6659	6.0379
	${2048}^2$	6.7595	46.7727	36.4367	29.9418	51.4417	29.3363	22.3794	28.6397	84.7245

TMAOR: two-step modulus-based accelerated over-relaxation; MSOR: modulus-based successive over-relaxation; MBIFP: modulus-based inexact fixed-point; CBB: Cauchy–Barzilai–Borwein.

Table 2.

The numerical results on Example 2.

$\mathbf{I}\mathbf{n}\mathbf{d}$	$\mathbf{n}$	MBIFP-CBB	MSOR ${}_a$	MSOR ${}_b$	MSOR ${}_c$	TMSOR ${}_a$	TMSOR ${}_b$	TMSOR ${}_c$	TMAOR ${}_a$	TMAOR ${}_b$
iter	${32}^2$	8	91	59	40	44	27	23	26	90
	${128}^2$	9	103	67	45	50	32	25	31	102
	${512}^2$	9	115	74	51	57	36	26	35	115
	${2048}^2$	10	127	82	56	63	40	29	40	127
avg-init	${32}^2$	1.1	1	1	1	2	2	2	2	2
	${128}^2$	1.2	1	1	1	2	2	2	2	2
	${512}^2$	1.3	1	1	1	2	2	2	2	2
	${2048}^2$	1.4	1	1	1	2	2	2	2	2
res	${32}^2$	7.7 $\times {10}^{-5}$	9.0 $\times {10}^{-5}$	9.5 $\times {10}^{-5}$	9.5 $\times {10}^{-5}$	8.2 $\times {10}^{-5}$	9.8 $\times {10}^{-5}$	8.6 $\times {10}^{-5}$	9.2 $\times {10}^{-5}$	9.7 $\times {10}^{-5}$
	${128}^2$	3.0 $\times {10}^{-5}$	9.1 $\times {10}^{-5}$	8.7 $\times {10}^{-5}$	9.8 $\times {10}^{-5}$	9.1 $\times {10}^{-5}$	7.5 $\times {10}^{-5}$	7.1 $\times {10}^{-5}$	8.0 $\times {10}^{-5}$	9.9 $\times {10}^{-5}$
	${512}^2$	8.9 $\times {10}^{-5}$	9.5 $\times {10}^{-5}$	9.9 $\times {10}^{-5}$	8.5 $\times {10}^{-5}$	8.1 $\times {10}^{-5}$	7.9 $\times {10}^{-5}$	9.7 $\times {10}^{-5}$	9.3 $\times {10}^{-5}$	9.2 $\times {10}^{-5}$
	${2048}^2$	2.0 $\times {10}^{-5}$	1.0 $\times {10}^{-4}$	9.6 $\times {10}^{-5}$	9.5 $\times {10}^{-5}$	8.9 $\times {10}^{-5}$	8.2 $\times {10}^{-5}$	7.0 $\times {10}^{-5}$	7.9 $\times {10}^{-5}$	9.6 $\times {10}^{-5}$
cput	${32}^2$	0.0010	0.0062	0.0039	0.0028	0.0056	0.0055	0.0033	0.0029	0.0110
	${128}^2$	0.0136	0.0970	0.0748	0.0495	0.1042	0.0765	0.0813	0.1036	0.2952
	${512}^2$	0.4322	3.8658	2.4643	1.7733	5.3265	2.1896	1.8659	2.6839	6.8261
	${2048}^2$	9.0873	61.8585	39.5533	26.9676	52.3150	34.0526	30.1440	41.3632	119.2776

MSOR: modulus-based successive over-relaxation; TMSOR: two-step modulus-based successive over-relaxation; TMAOR: two-step modulus-based accelerated over-relaxation;MBIFP: modulus-based inexact fixed-point; CBB: Cauchy–Barzilai–Borwein.

Example 1

In SNCP (1), $A$ is symmetric and $f(u)=Au+\varphi (u)+q$ , where $A=\hat{A}+\mu I\in {\mathbb{R}}^{n\times n}$ ,

$$\begin{array}{rl}& \hat{A}=\left(\begin{array}{cccc}H& -I& & \\ -I& H& \ddots & \\ & \ddots & \ddots & -I\\ & & -I& H\end{array}\right)\text{with}\\ & H=\left(\begin{array}{cccc}4& -1& & \\ -1& 4& \ddots & \\ & \ddots & \ddots & -1\\ & & -1& 4\end{array}\right)\mathrm{and}\mu =2,\end{array}$$

and $\varphi (u)=(u_1/(1+u_1),\dots ,u_n/(1+u_n){)}^T$ , $q=(-1,1,\dots ,(-1{)}^n{)}^T$ .

Example 2

In SNCP (1), $A$ is asymmetric and $f(u)=Au+\varphi (u)+q$ , where $A=\hat{A}+\mu I\in {\mathbb{R}}^{n\times n}$ ,

$$\begin{array}{rl}& \hat{A}=\left(\begin{array}{cccc}H& -0.5I& & \\ -1.5I& H& \ddots & \\ & \ddots & \ddots & -0.5I\\ & & -1.5I& H\end{array}\right)\text{with}\\ & H=\left(\begin{array}{cccc}4& -0.5& & \\ -1.5& 4& \ddots & \\ & \ddots & \ddots & -0.5\\ & & -1.5& 4\end{array}\right)\mathrm{and}\mu =4,\end{array}$$

and $\varphi (u)=(\arctan (u_1),\dots ,\arctan (u_n){)}^T$ , $q=(1,-1,\dots ,(-1{)}^n{)}^T$ .

Remark 2

(i) The nonlinear terms $\varphi (\cdot )$ in Examples 1-2 are ${\ell }_{\varphi }$ -Lipschitz with ${\ell }_{\varphi }=1$ . (ii) Figure 2 validates the reasonability of the second claim of condition (2) with respect to Example 1 and Example 2, where the dimension of variable $u$ varies from ${10}^2$ to ${100}^2$ .

Figure 2.

The relationship between $\delta ={\ell }_{\varphi }+{\lambda }_{A^N}^{max}$ and ${\lambda }_1^S$ associated to Examples 1 and 2.

The proposed MBIFP-CBB method is compared with MSOR, TMSOR and TMAOR methods on Examples 1 and 2. The method that has the least CPU time (in seconds) is selected as the winner and represented by bold fonts. Besides, to further validate the performance of these methods, we set $\varphi (u)=u/(1+u)$ , $\arctan (u)$ , $\ln (1+u)$ and $(u-\sin (u))/2$ , respectively. The MBIFP-CBB method is also compared to MSOR ${}_b$ , MSOR ${}_c$ , TMSOR ${}_b$ and TMAOR ${}_a$ methods when they are used to Examples 1 and 2 by setting $x_0=ones(n,1)$ . The results that the logarithm of residuals vs CPU time (in ${10}^{-2}$ seconds) are displayed in Figures 3 and 4.

Figure 3.

Residual comparision with $x_0=ones(n,1)$ for Example 1 (symmetric structured nonlinear complementarity problem (SNCP)).

Figure 4.

Residual comparision with $x_0=ones(n,1)$ for Example 2 (non-symmetric structured nonlinear complementarity problem (SNCP)).

From Tables 1 and 2 and Figures 3 and 4, one can find that all methods are convergent in the numerical experiments. In almost all settings in these experiments, the number of total iterations and the CPU time of the MBIFP-CBB method have the most obvious superiority to that of the MSOR, TMSOR, and TMAOR methods.

Conclusions

The SNCP is reformulated as a nonsmooth implicit fixed-point system $x=H(x)$ founded on the modulus decomposition. Then, we have been developed a MBIFP iteration method for the solution of SNCP. The global convergence of the proposed method mainly relies on a crucial condition that $H(x)$ is a contraction mapping. Hence, a sufficient condition on parameter $\omega$ is analyzed for the contractibility of $H(x)$ . Some preliminary numerical results show that, compared with some existing matrix splitting type methods, the proposed MBIFP method cooperated with the CBB method is applicable and efficient for large-scale SNCP. It is worth noting that the proposed MBIFP method provides a general framework for solving the SNCP. Therefore, we can select an effective linear solver for the iteration subproblem by utilizing the structure of the system matrix $A$ in SNCP.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received the following financial support for the research, authorship, and/or publication of this article: This work was supported partly by the National Natural Science Foundation of China with grant 12071398, the Natural Science Foundation of Hunan Province with grant 2020JJ4567 and the Key Scientific Research Found of Hunan Education Department with grants 20A097.

ORCID iD

Zheng Peng