A smooth Levenberg-Marquardt method without nonsingularity condition for wLCP

1.   Introduction

This paper considers the weighted Linear Complementarity Problem (wLCP) introduced by Potra ^[14] which consists in finding vectors $x\in\mathbb{R}^n, s\in\mathbb{R}^n, y\in\mathbb{R}^m$ such that

Here $P\in \mathbb{R}^{(n+m)\times n}, Q\in \mathbb{R}^{(n+m)\times n}, R\in \mathbb{R}^{(n+m)\times m}$ are given matrices, $d\in \mathbb{R}^{n+m}$ is a given vector, $w\geq0$ is a given weight vector (the data of the problem) and ${xs}$ is the componentwise product of the vectors $x$ and $s$. The matrix $R$ is assumed to have full column rank. If the weight vector $w$ is chosen to be the zero vector, then the wLCP reduces to the general LCP studied in ^[21]. The wLCP is called monotone, if

The wLCP can be used for modeling a larger class of problems from science and engineering. For example, Fisher's competitive market equilibrium model can be formulated as a wLCP that can be efficiently solved by interior-point methods ^[14]. Moreover, the Quadratic Programming and Weighted Centering problem, which generalizes the notion of a Linear Programming and Weighted Centering problem proposed by Anstreicher ^[1], can be formulated as a monotone wLCP ^[14]. Lately, Chi et al. ^[3] and Gowda ^[10] studied wLCP on Euclidean Jordan algebras.

Since Potra introduced the notion of wLCP ^[14], many numerical algorithms have been proposed for solving this problem. One class of effective algorithms is interior-point methods. For example, Potra ^[14] proposed two interior-point methods for solving general monotone wLCPs. In ^[15], Potra proposed a corrector-predictor interior-point method for solving the sufficient wLCP. Asadi et al. ^[2] introduced a full-Newton step interior-point method for solving the monotone wLCP. Chi et al. ^[4] proposed a full-modified-Newton infeasible interior-point method for solving a special wLCP. It is worth pointing out that, to establish the computational complexity of interior-point methods, one usually requires that the wLCP is monotone (e.g., ^[2,14]). Another class of effective algorithms for solving the wLCP is smoothing Newton-type algorithms. This class of algorithms is to use a smoothing function to reformulate the wLCP as a system of smooth nonlinear equations and then solve it by Newton method. For example, Zhang ^[22] presented a smoothing Newton algorithm for solving the wLCP. Tang ^[16] proposed a variant nonmonotone smoothing algorithm for solving the wLCP with improved numerical results. Tang and Zhang ^[17] proposed a nonmonotone smoothing Newton algorithm for solving general wCPs. Notice that, to ensure Newton step be feasible, smoothing Newton-type algorithms in ^[16,17,22] also require that the wLCP is monotone. Moreover, to obtain local fast convergence rate, smoothing Newton-type algorithms in ^[16,22] need the nonsingularity condition.

Lately, based on a nonsmooth weighted complementarity function, Tang and Zhou ^[19] reformulated the wLCP as a nonsmooth nonlinear equation and proposed a damped Gauss-Newton method to solve it. Their method can be used to solve nonmonotone wLCPs and has local quadratic convergence under the local error bound condition which is weaker than the nonsingularity condition. Motivated by their work, in this paper we introduce a weighted complementarity function which is smooth everywhere. By using this function, we reformulate the wLCP as a smooth nonlinear equation and propose a Levenberg-Marquardt method to solve it. Different from current Levenberg-Marquardt type methods (e.g., ^[7,8,9]), the proposed method is designed based on a simple derivative-free line search technique. Compared with smoothing Newton-type algorithms in ^[16,17,22], the proposed method has two advantages. (i) It is well-defined and is globally convergent without any additional condition. Hence it can be used to solve nonmonotone wLCPs. (ii) It has local sub-quadratic convergence rate under the local error bound condition. It is worth pointing out that, to obtain the local fast convergence rate, classical Levenberg-Marquardt methods (e.g., ^[7,8,9]) also need the condition that the Jacobian is Lipschitz continuous. In this paper we show that this condition holds for our method (see, Lemma 4.1 below). We also report some numerical results which indicate that our method is more effective for solving monotone and nonmonotone wLCPs than the damped Gauss-Newton method studied in ^[19].

This paper is organized as follows. In Section 2, we reformulate the wLCP as a smooth nonlinear equation and propose a Levenberg-Marquardt method to solve it. In Section 3, we give its global convergence. In Section 4, we analyze its local sub-quadratic convergence under the local error bound condition. Numerical results are reported in Section 5. Some conclusions are given in Section 6.

Throughout this paper, $\mathbb{R}^n$ denotes the set of all $n$ dimensional real vectors. All vectors are column vectors and for simplicity, the column vector $(u_1^T, \ldots, u_n^T)^T$ is written as $(u_1, \ldots, u_n)$ where $u_i\in\mathbb{R}^{n_i}$. $\|\cdot\|$ denotes the 2-norm. For any vector $x\in {\mathbb{R}}^n$, we denote $x$ by $\mathrm{vec}(x_i)$ and the diagonal matrix whose $i$th diagonal element is $x_i$ by $\mbox{diag}(x_i)$. For a given set $S\subset \mathbb{R}^n$ and for any $u\in \mathbb{R}^n$, $\mbox{dist}(u, S) = \inf\limits_{v\in S}\{\|u-v\|\}$. For any $\alpha, \beta > 0,$ $\alpha = {O} (\beta)$ (respectively, $\alpha = o(\beta))$ means that $\lim\sup_{\beta\rightarrow0}\frac{\alpha}{\beta} < \infty$ (respectively, $\lim\sup_{\beta\rightarrow0}\frac{\alpha}{\beta} = 0$).

2.   A smooth Levenberg-Marquardt method

2.1.   The reformulation of wLCP

The weighted complementarity function serves an important role in designing Newton-type methods for solving the wLCP. For a fixed $c\geq0,$ a function $\phi^c(a, b): \mathbb{R}^2\rightarrow \mathbb{R}$ is called a weighted complementarity function if it satisfies

One popular weighted complementarity function is

When $c = 0,$ $\phi^c(a, b)$ is the well-known Fischer-Burmeister function for nonlinear complementarity problems. Obviously, $\phi^c(a, b)$ is not continuously differentiable (smooth) everywhere. By using $\phi^c(a, b)$, Tang and Zhou ^[19] reformulated the wLCP as the following nonsmooth nonlinear equation:

where $w = (w_1, ..., w_n)^T\geq0$ is the weight vector given in the wLCP. Tang and Zhou ^[19] presented a damped Gauss-Newton method to solve $\Phi (x, s, y) = 0$ and established its local quadratic convergence under the local error bound condition.

In this paper, for a fixed $c\geq0,$ we consider the following nonnegative function:

The following lemma shows that $\psi^c$ is a weighted complementarity function and it is smooth everywhere.

Lemma 2.1. Let $\psi^c$ be defined by $(2.1)$.Then the following results hold.

(i) $\psi^c$ is a weighted complementarity function, i.e.,

(ii) $\psi^c$ is continuously differentiable at any $(a, b)\in \mathbb{R}^2$ and

where $\nabla_a\psi^c(0, 0) = \nabla_b\psi^c(0, 0) = -\sqrt{2c}$ and for any $(a, b)\neq(0, 0),$

(iii) For any $(a, b)\in \mathbb{R}^2$, one has

(iv) For any $(a, b)\in \mathbb{R}^2$, one has

Proof. The result (i) obviously holds. By a direct computation, we can obtain the result (ii). Moreover, for any $(a, b)\in \mathbb{R}^2$, by the result (ii), we have $\nabla_a\psi^c(0, 0)\nabla_b\psi^c(0, 0) = 2c\geq0$, and if $(a, b)\neq(0, 0)$, then we also have

where the inequality holds because

Now we prove the result (iv). First, we prove the implication

Since $“\Longrightarrow”$ obviously holds, we only need to prove $“\Longleftarrow”$. By the result (ii), $\nabla_a\psi^c(0, 0) = 0$ gives $c = 0$ and so $\psi^c(0, 0) = 0$. For $(a, b)\neq(0, 0)$, $\nabla_a\psi^c(a, b) = 0$ gives $\phi^c(a, b) = 0$ or $1-\frac{a}{\sqrt{a^2+b^2+2c}} = 0.$ If $\phi^c(a, b) = 0$, then $\psi^c(a, b) = 0$ and we are done. If $1-\frac{a}{\sqrt{a^2+b^2+2c}} = 0,$ then $a = \sqrt{a^2+b^2+2c}$ which yields $b = c = 0$ and $a = |a|$. These implies that $\phi^c(a, b) = 0$ and so $\psi^c(a, b) = 0$. By the same way, we can also prove $\psi^c(a, b) = 0\Longleftrightarrow\nabla_b\psi^c(a, b) = 0.$ The proof is completed.$

Figures 1 and 2 illustrate the geometrical interpretations of $\phi^c$ and $\psi^c$ with $c = 1$ which show that they are very different.

By using the weighted complementarity function $\psi^c$ given in (2.1), we now reformulate the wLCP as the following smooth nonlinear equation:

By Lemma 2.1, we have the following result.

Lemma 2.2. (i) $\mathrm{H} (x, s, y) = 0$ if and only if $(x, s, y)$ is a solution of the ${\rm{wLCP}}$.

(ii) $\mathrm{H} (x, s, y)$ is continuously differentiable at any $(x, s, y)\in\mathbb{R}^{2n+m}$ withthe Jacobian

in which $\nabla_{x_i}\psi^{w_i}(\cdot, \cdot)$ and $\nabla_{s_i}\psi^{w_i}(\cdot, \cdot)$ are given in Lemma 2.1 (ii).

2.2.   The algorithm

Let $z: = (x, s, y)$ and $\mathrm{H} (z)$ be given in (2.2). We define the merit function $\Psi:\mathbb{R}^{2n+m}\to \mathbb{R}$ as

Then, by Lemma 2.2, we have the following result.

Lemma 2.3. (i) $\Psi(z) = 0$ if and only if $z = (x, s, y)$ is a solution of the ${\rm{wLCP}}$.

(ii) $\Psi(z)$ is continuously differentiable at any $z\in\mathbb{R}^{2n+m}$ with its gradient $\nabla \Psi(z) = \mathrm{H} '(z)^T\mathrm{H} (z)$.

Now we describe our method as follows.

Algorithm 2.1. (A smooth Levenberg-Marquardt method)

Step 0: Choose parameters $\theta, \rho, \gamma\in(0, 1)$ and an initial point $z^0: = (x^0, s^0, y^0)$.Set $k: = 0.$

Step 1: If $\nabla \Psi(z^k) = 0$, then stop.

Step 2: Set $\mu_k: = \theta\|\mathrm{H}(z^k)\|^\delta$ where $\delta\in [1, 2]$ is a constant. Compute the search direction $d_k\in \mathbb{R}^{2n+m}$ by solving

Step 3: Find a step-size $\alpha_k : = \rho^{m_k}$, where $m_k$ is the smallest nonnegative integer $m$ satisfying

Step 4: Set $z^{k+1}: = z^k+\alpha_kd_k$. Set $k: = k+1$ and go to Step 1.

The algorithmic framework of Algorithm 2.1 is much simpler than many Levenberg-Marquardt type methods (e.g, ^{[7,8,9,13,18,23]}). The main feature of Algorithm 2.1 is that it adopts a simple derivative-free line search in Step 3 which avoids computing the gradient $\nabla \Psi(z^k)$.

Theorem 2.1. Algorithm $2.1$ is well-defined.

Proof. For some $k$, if $\nabla \Psi(z^k)\neq0$, then $\mathrm{H}(z^k)\neq0$ and hence $\mu_k = \theta\|\mathrm{H}(z^k)\|^\delta > 0$. So, $\mathrm{H}'(z^k)^T\mathrm{H}'(z^k)+\mu_kI$ is positive definite and the search direction $d_k$ in Step 2 is well-defined. Since $\nabla \Psi(z^k)\neq0$, by (2.5) we have $d_k\neq0$ and

This implies that $d_k$ is a descent direction of the merit function $\Psi(z)$ at $z^k.$ Next, we show that there exists at least a nonnegative integer $m$ satisfying (2.6). On the contrary, we suppose that for any nonnegative integer $m$,

Multiplying both sides of the above inequality by $\frac{1}{2}[\|\mathrm{H}(z^k+\rho^{m} d_k)\|+\|\mathrm{H}(z^k)\|]$, we have

Since $\Psi$ is continuously differentiable at $z^k$, by letting $m\rightarrow \infty$ in (2.8), we have $\nabla\Psi(z^k)^Td_k\geq0$ which contradicts (2.7). So, we can find a step-size $\alpha_k$ in Step 3 and get the $(k+1)$-th iteration $z^{k+1} = z^k+\alpha_kd_k$ in Step 4. Hence, Algorithm 2.1 is well-defined.

3.   Global convergence

In the following, we assume $\nabla \Psi(z^k)\neq0$ for all $k\geq0$, so that Algorithm 2.1 generates an infinite sequence $\{z^k\}$. To establish the global convergence, we need the following lemma.

Lemma 3.1. Let $\{z^k\}$ be the sequence generated by Algorithm $2.1$.Then one has:

(i) $\|d_k\|\leq\frac{1}{2\sqrt{\theta}}\|\mathrm{H}(z^k)\|^{1-\frac{\delta}{2}}$ for all $k\geq0;$

(ii) $\lim\limits_{k\rightarrow \infty}\|\alpha_kd_k\| = 0.$

Proof. For any $k\geq0,$ we suppose that the singular value decomposition of $\mathrm{H}'(z^k)$ is

where $U_k$ and $V_k$ are orthogonal matrices and $\Sigma_k = {{\mathrm{diag}}}(\sigma_1^k, ..., \sigma_r^k, 0, ..., 0)$ with $\sigma_1^k\geq \cdots\geq\sigma_r^k > 0.$ Then, by (2.5) we have for any $k\geq0$,

where the inequality holds because $\frac{\sigma_i^k}{(\sigma_i^k)^2+\mu_k}\leq\frac{1}{2\sqrt{\mu_k}}$ for all $i = 1, ..., r$. This together with $\mu_k = \theta\|\mathrm{H}(z^k)\|^\delta$ gives the result (i). Moreover, by (2.6) we have

Since $\{\|\mathrm{H}(z^k)\|\}$ is a monotonically decreasing sequence by (2.6), there exists a constant $\mathrm{H}^*\geq 0$ such that $\lim\limits_{k\rightarrow \infty}\|\mathrm{H}(z^k)\| = \mathrm{H}^*$. This together with (3.1) proves the result (ii).

Theorem 3.1. Let $z^*$be any accumulation point of the sequence $\{z^k\}$ generated by Algorithm $2.1$. Then$z^*$ is a stationary point of the merit function $\Psi(z)$, i.e., $\nabla\Psi(z^*) = 0$.Moreover, if $\mathrm{H}'(z^*)$ is nonsingular, then $\mathrm{H}(z^*) = 0$ and so $z^* = (x^*, s^*, y^*)$ is a solution of the ${rm{wLCP}}$.

Proof. Without loss of generality, we may assume that $z^*$ is the limit of the subsequence $\{z^k\}_{k\in K}$ where $K\subset\{0, 1, ...\}$, i.e, $\lim\limits_{(K\ni)k\rightarrow \infty}z^k = z^*.$ Then, by the continuity of $\mathrm{H}$ and $\mathrm{H}'$,

and consequently

Obviously, $\nabla\Psi(z^*) = 0$ when $\mathrm{H}(z^*) = 0$. Now we assume $\|\mathrm{H}(z^*)\| > 0$ and will derive a contradiction. Since $\{\|\mathrm{H}(z^k)\|\}$ is a monotonically decreasing sequence, by Lemma 3.1 (i), we have for all $k\in K$,

Thus, $\{d_k\}_{k\in K}$ has a convergent subsequence and we may assume $\lim\limits_{(K_1\ni)k\rightarrow \infty}d_k = d^*$ where $K_1\subset K.$ In the following, we show that $d^* = 0$. In fact, if $d^*\neq0$, then by Lemma 3.1 (ii), we have $\lim\limits_{(K_1\ni)k\rightarrow \infty}\alpha_k = 0$. Moreover, from Step 3 we have for all $k\in K_1,$

Multiplying both sides of this inequality by $\frac{1}{2}\big[\|\mathrm{H}(z^k+\rho^{-1}\alpha_kd_k)\|+\|\mathrm{H}(z^k)\|\big]$, we have for all $k\in K_1,$

Since $\Psi$ is continuously differentiable at $z^*$, by letting $k\rightarrow \infty$ with $k\in K_1$ in (3.2), we have

On the other hand, by (2.7) we have

So, from (3.3) and (3.4) we have $\nabla\Psi(z^*)^Td^* = 0$, which together with (2.5) gives

Since $\|\mathrm{H}(z^*)\| > 0,$ the matrix $\mathrm{H}'(z^*)^T\mathrm{H}'(z^*)+\theta\|\mathrm{H}(z^*)\|^\delta I$ is positive definite. Hence, we have $d^* = 0$. A contradiction is derived. Therefore, we have $d^* = 0.$ Then, by letting $k\to \infty$ with $k\in K_1$ in (2.5), we have

This proves the first result. The second result follows from $\nabla \Psi(z^*) = \mathrm{H} '(z^*)^T\mathrm{H} (z^*) = 0$. We complete the proof.

Theorem 3.2. Let $\{z^k\}$ be the sequence generated by Algorithm $2.1$. If $\{z^k\}$ has an isolated accumulation point$z^*$, then the whole sequence $\{z^k\}$ converges to $z^*$.

Proof. By Lemma 3.1 (ii), we have $\lim\limits_{ k\rightarrow \infty}\|z^{k+1}-z^k\| = 0.$ This together with ^[6,Proposition 8.3.10] proves the theorem.

At the end of this section, we consider a special wLCP which consists in finding vectors $x\in\mathbb{R}^n, s\in\mathbb{R}^n$ such that

where $q\in \mathbb{R}^n$ and $M\in \mathbb{R}^{n\times n}$ is a $P_0$ matrix, that is, for any $\xi\in \mathbb{R}^n$ with $\xi\neq0$, there exists an index $i_0\in \{1, ..., n\}$ such that $\xi_{i_0}\neq0$ and $\xi_{i_0}(M\xi)_{i_0}\geq0$.

When $w$ is chosen to be the zero vector, the $P_0$ wLCP (3.5) reduces to the following standard $P_0$ LCP:

which has many applications in economics and engineering and has been extensively studied in literatures (e.g., ^{[11,12,24,25,26]}). By using the weighted complementarity function $\psi^c(a, b)$ defined by (2.1), we can reformulate the $P_0$ wLCP as the smooth nonlinear equation

and apply Algorithm 2.1 to solve it. For the $P_0$ wLCP, we have the following global convergence result.

Theorem 3.3. Let $\{(x^k, s^k)\}$ be the sequence generated by Algorithm $2.1$ for solving the nonlinear equation $(3.6)$.Then any accumulation point $(x^*, s^*)$ of $\{(x^k, s^k)\}$ is a solution of the $P_0$ ${\rm{wLCP}}$.

Proof. Since the $P_0$ wLCP is a special case of the wLCP, Theorem 3.1 still holds. So, we have $\nabla \Psi(x^*, s^*) = \mathrm{H} '(x^*, s^*)^T\mathrm{H} (x^*, s^*) = 0$. By (3.6), we have

in which $\nabla_{x_i}\psi^{w_i}(\cdot, \cdot)$ and $\nabla_{s_i}\psi^{w_i}(\cdot, \cdot)$ are given in Lemma 2.1 (ii). Then, it follows from $\mathrm{H} '(x^*, s^*)^T\mathrm{H} (x^*, s^*) = 0$ that

Now we assume $Mx^*+q-s^*\neq0$. Since $M$ is a $P_0$ matrix and so is $M^T,$ there exists an index $i_0\in \{1, ..., n\}$ such that $(Mx^*+q-s^*)_{i_0}\neq0$ and

By (3.8) and (3.9), we have

which together with (3.10) gives

On the other hand, by Lemma 2.1 (iii), we have

From (3.13) and (3.14) it holds that

This together with Lemma 2.1 (iv) gives $\psi^{w_{i_0}}(x^*_{i_0}, s^*_{i_0}) = 0$. Then, by (3.12) we have $(Mx^*+q-s^*)_{i_0} = 0$ which contradicts the choice of the index $i_0$ such that $(Mx^*+q-s^*)_{i_0}\neq0$. Hence, $Mx^*+q-s^* = 0$. Furthermore, by (3.8) we have $\nabla_{x_i}\psi^{w_i}(x^*_i, s^*_i)\psi^{w_i}(x^*_i, s^*_i) = 0$ which together with Lemma 2.1 (iv) yields $\psi^{w_i}(x^*_i, s^*_i) = 0$ for all $i = 1, ..., n$. Therefore, $\mathrm{H}(x^*, s^*) = 0$ and $(x^*, s^*)$ is a solution of the $P_0$ wLCP.

4.   Local sub-quadratic convergence

Let $\mathrm{H}(z)$ be given in (2.2). Denote

In this section, we assume that the whole sequence $\{z^k\}$ generated by Algorithm 2.1 converges to some point $z^*\in\mathcal{Z}^*$. Then $\lim\limits_{k\rightarrow \infty}\mathrm{H}(z^k) = \mathrm{H}(z^*) = 0.$ Now we make the following assumption.

Assumption 4.1. $\mathrm{H}(z)$ provides a local error bound on some neighbourhood of $z^*$, i.e., thereexist constants $\xi > 0$ and $\epsilon > 0$ such that

As it is well-known, Assumption 4.1 is the local error bound condition which is weaker than the nonsingularity condition. It is worth pointing out that, to obtain the local quadratic convergence, classical Levenberg-Marquardt methods (e.g., ^[7,8,9]) also need to assume that the Jacobian is Lipschitz continuous on the set $N(z^*, \epsilon)$. In the following, we show that this assumption holds for our method.

Lemma 4.1. The Jacobian $\mathrm{H}'(z)$ given in $(2.3)$ is Lipschitz continuous on $\mathbb{R}^{2n+m}$, i.e., there exists a constant $M > 0$ such that

Proof. Obviously, we only need to prove that the gradient $\nabla\psi^c(a, b)$ given in Lemma 2.1 (ii) is Lipschitz continuous on $\mathbb{R}^2$. First, by ^[5,Lemma 3.1], $\nabla\psi^c(a, b)$ is Lipschitz continuous on $\mathbb{R}^2$ when $c = 0.$ Now we consider $c > 0$. Then $\psi^c(a, b)$ is twice continuously differentiable at any $(a, b)\in \mathbb{R}^2$ with

where

Since

also notice that

we have

Moreover, it is clear that

Thus, we have

This implies that there exists a constant $C > 0$ independent of $(a, b)$ such that

Then, by the Mean-Value Theorem, we have

holds for all $(a, b), (\tilde{a}, \tilde{b})\in \mathbb{R}^2$. The proof is completed.

By Lemma 4.1, we can directly have

and there exists a constant $L > 0$ such that

In the following, we denote $\bar{z}^k$ as the vector in $\mathcal{Z}^*$ that satisfies

Lemma 4.2. Let $\{z^k\}$ be the sequence generated by Algorithm $2.1$.If Assumption $4.1$ holds, then for all sufficiently large $k,$

Proof. Notice that for all sufficiently large $k,$

which implies that $\bar{z}^k$ sufficiently close to $z^*$. So, by (4.5), for all sufficiently large $k,$

For any $k\geq0$, we consider the following optimization problem:

Then, the search direction $d_k$ generated by (2.5) is the solution of (4.9) because $\varphi_k(d)$ is a strictly convex quadratic function and $d_k$ is a stationary point of $\varphi_k(d)$. Hence, for all sufficiently large $k,$ by (4.4) and (4.8) we have

This together with (4.9) yields

Moreover, for all sufficiently large $k,$ by Assumption 4.1, we have

which together with (4.9) and (4.10) gives

The proof is completed.

Theorem 4.1. Let $\{z^k\}$ be the sequence generated by Algorithm $2.1$.If Assumption $4.1$ holds, then for all sufficiently large $k,$ we have

Proof. By (4.7), for all sufficiently large $k,$

where $c_1: = 1+\sqrt{\frac{M^2+\theta L^\delta}{\theta\xi^\delta}}$. This implies that $z^k+d_k$ sufficiently close to $z^*$. Hence, by (4.4), for all sufficiently large $k,$

which together with (4.6) and (4.7) gives

where $c_2: = \sqrt{M^2+\theta L^\delta}+{\frac{M(M^2+\theta L^\delta)}{\theta\xi^\delta}}.$ Thus, by (4.1) and (4.15), for all sufficiently large $k,$

Let ${\tilde{z}^k}$ be the vector in $\mathcal{Z}^*$ that satisfies $\|{{z^k+d_k}}-{\tilde{z}^k}\| = {\mathrm{dist}}({{z^k+d_k}}, \mathcal{Z}^*).$ Then, by (4.13) and (4.16), for all sufficiently large $k,$

which yields $\tilde{z}^k$ sufficiently close to $z^*$. So, by (4.1), (4.5) and (4.16), also notice that $\mathrm{H}({\tilde{z}^k}) = 0$ as ${\tilde{z}^k}\in\mathcal{Z}^*$, for all sufficiently large $k,$

Moreover, by (4.1) and (4.7), for all sufficiently large $k,$

Hence, by (4.17) and (4.18) we have

which implies that for all sufficiently large $k,$

This shows that the step-size $\alpha_k = 1$ is accepted in Step 3 for all sufficiently large $k$. Consequently, for all sufficiently large $k,$ we have $z^{k+1} = z^k+d_k$ which together with (4.16) prove the theorem.

By Theorem 4.1, similarly as the proof of ^[19,Theorem 5.2], we can obtain the following sub-quadratic convergence property.

Theorem 4.2. Let $\{z^k\}$ be the sequence generated by Algorithm $2.1$.If Assumption $4.1$ holds, then one has

5.   Numerical results

In this section, we report some numerical results of Algorithm 2.1. All experiments are carried on a PC with CPU of Inter(R) Core(TM)i7-7700 CPU @ 3.60 GHz and RAM of 8.00GB. The codes are written in MATLAB and run in MATLAB R2018a environment. The parameters used in Algorithm 2.1 are chosen as $\theta = 10^{-4}, \rho = 0.8, \gamma = 10^{-4}, \delta = 1.$

We apply Algorithm 2.1 to solve the wLCP (1.1) in which

where $A \in \mathbb{R}^{m \times n}$ is a full row rank matrix with $m < n$, $M\in \mathbb{R}^{n \times n}$, $b\in\mathbb{R}^m$ and $f\in\mathbb{R}^n.$ This example comes from ^[14]. For the purposes of comparison, we also apply the the damped Gauss-Newton method studied by Tang and Zhou ^[19] to solve this test problem. In our experiments, we test the following two class of wLCPs.

(The monotone wLCP) We choose $A = \bf{randn} (m, n)$ with the rank of $A$ being $m$ and $M = BB^T/\|BB^T\|$ with $B = \bf{rand}(n, n)$. Then we choose $\hat{x} = \bf{rand}(n, 1)$, $f = \bf{rand}(n, 1)$ and set $b: = A\hat{x}$, $\hat{s}: = M\hat{x}+f$ and $w: = \hat{x}\hat{s}$. This wLCP is monotone. For each problem with sizes $n( = 2m)$, we generate ten instances and solve them by using the following three starting points:

(i) $x^0 = s^0 = (1, ..., 1)^T$, $y^0 = (0, ..., 0)^T$;

(ii) $x^0 = s^0 = (1, 0, ..., 0)^T$, $y^0 = (0, ..., 0)^T$; (iii) $x^0 = \bf{rand}(n, 1), s^0 = \bf{rand}(n, 1), y^0 = \bf{rand}(m, 1)$.

We use $\|{\mathrm{H}}(z^k)\|\leq10^{-5}$ as the stopping criterion. Numerical results are listed in Table 1 where SLMM denotes the smooth Levenberg-Marquardt method studied in this paper, DGNM denotes the damped Gauss-Newton method studied in ^[19], SP denotes the starting point, AIT and ACPU denote the average number of iterations and the average CPU time in seconds respectively. From Table 1, we can see that SLMM has the advantage over DGNM, especially for large scale test problems.

(The nonmonotone wLCP) We choose $M = B1/\|B1\|-B2/\|B2\|$ with $B1 = \bf{rand}(n, n)$ and $B2 = \bf{rand}(n, n)$. The matrix $A$ and vectors $b, f, w$ are generated by the same way as before. Since the matrix $M$ is not symmetric positive semidefinite, this class of wLCP may be nonmonotone. For each problem with sizes $n( = 2m)$, we also generate ten instances and solve them by using three starting points as before. Moreover, we use $\|{\mathrm{H}}(z^k)\|\leq10^{-5}$ and $iter < 50$ as the stopping criterion where iter denotes the number of iterations. Numerical results are listed in Table 2, where $*$ stands for that the method fails to solve some instances as the iteration number is grater than 50 and the average is based on the successful instances through our numerical reports. From Table 2, we can see that, although both SLMM and DGNM can be applied to solve nonmonotone wLCPs, the former has better numerical performance than the latter.

6.   Conclusions

Based on a smooth weighted complementarity function, we reformulated the wLCP as a smooth nonlinear equation and proposed a Levenberg-Marquardt method to solve it. The proposed method is well-defined and it is globally convergent without any additional condition. Moreover, we proved that the proposed method has local sub-quadratic convergence rate under the local error bound condition which is weaker than the nonsingularity condition. Numerical results show that our method is very effective for solving monotone and nonmonotone wLCPs. In Algorithm 2.1, the Eq (2.5) is solved exactly in each iteration which maybe expensive for large-scale wLCPs. As a future research issue, it is worth investigating smooth inexact Levenberg-Marquardt method for solving wLCPs. Moreover, Wang and Fan ^[20] lately established the local convergence rate of Levenberg-Marquardt method under the Hölderian local error bound condition which is more general than Assumption 4.1 in this paper. Thus, another interesting issue is whether Algorithm 2.1 in this paper has local fast convergence rate under the Hölderian local error bound condition.

Acknowledgments

This paper is partly supported by Natural Science Foundation of Henan Province (222300420520) and Key scientific research projects of Higher Education of Henan Province (22A110020).

Conflict of interest

All authors declare no conflicts of interest in this paper.