General fixed-point method for solving the linear complementarity problem

1.   Introduction

In this paper, we consider the linear complementarity problem, which is to find a vector $x\in {\mathbb{R}}^{n}$ such that

where $A = (a_{ij})\in{\mathbb{R}}^{n \times n}$ and $q\in{\mathbb{R}}^{n}$ are given. For convenience, such problem is usually abbreviated as LCP($A, q$), which is related to many practical problems, such as American option pricing problems, the market equilibrium problems and the free boundary problems for journal bearings, see ^[1,2,3,4,5] and the references therein.

To obtain the numerical solution of the LCP($A, q$) with a large and sparse system matrix, many kinds of iteration methods have been proposed and analyzed in recent decades. The projected method is a well known iteration method, which includes the smoothing projected gradient method ^[6], the projected gradient method ^[7], the partial projected Newton method ^[8], the improved projected successive over relaxation (IPSOR) ^[9], and so on. For more detailed materials about projected type iteration methods, we refer readers to ^{[10,11,12,13]} and the references therein. There are other efficient iteration methods for solving the LCP, e.g., the modulus-based matrix splitting iteration methods ^[2] and the nonstationary extrapolated modulus algorithms ^[14], the two-sweep modulus-based matrix splitting iteration methods ^[15], the general modulus-based matrix splitting method ^[5], the two-step modulus-based matrix splitting iteration method ^[16], the accelerated modulus-based matrix splitting iteration methods ^[17]. The main difference between the projected type iteration methods and the modulus-based type matrix splitting iteration methods is that the projected type iteration methods directly construct the iterative forms on an equivalent fixed-point equation based on the projection approaches and matrix splittings, however, the modulus-based type matrix splitting iteration methods reformulate the LCP($A, q$) as an implicit fixed-point equation by introducing positive diagonal parameter matrices, then construct the iterative forms based on all kinds of matrix splittings and avoid the projections. The most prominent difference is that one requires the projection and the other does not. For other iteration methods for solving the complementarity problems, we refer readers to ^{[18,19,20,21,22,23,24]} and the references therein.

In ^[4], Shi, Yang, and Huang presented a fixed-point (FP) method for solving a concrete LCP($A, q$) arising in American option pricing problems. The FP method is based on a fixed-point equation and belongs to the projected type iteration methods. However, the numerical examples in ^[4] shows that the number of iteration steps is very large when a suitable approximate solution is obtained. In this paper, we further discuss the projected type iteration methods and consider a general fixed-point equation by introducing a positive diagonal parameter matrix $\Omega$. We note that ^[4] considered the case where $\Omega = \alpha I$ with $\alpha > 0$ and $I$ is the identity matrix. We prove the equivalence between the general fixed-point equation and the linear complementarity problem. Based on the new fixed-point equation, we propose the general fixed-point (GFP) method with two iteration forms. We discuss the convergence conditions and provide the concrete convergence domains for the proposed method. Moreover, we discuss the optimal parameter problem and obtain an optimal parameter value.

The rest of this paper is organized as follows. In Section 2, we introduce some notations and concepts briefly, and then provide two lemmas which are required to derive the new fixed-point method. In Section 3, we propose the general fixed-point method with convergence analysis. In Section 4, we present numerical examples to illustrate the efficiency of the proposed method. Finally, we give the concluding remark in Section 5.

2.   Preliminaries

In this section, we first briefly review some notations and concepts, then provide two lemmas that will be used in Section 3.

A matrix $A = (a_{ij}) \in \mathbb{R}^{n\times n}$ is denoted by $A\geq 0\; (\mathrm{or}\; A > 0)$ if $a_{ij}\geq 0\; (\mathrm{or}\; a_{ij} > 0)$, and the absolute value matrix of $A = (a_{ij}) \in \mathbb{R}^{n\times n}$ is denoted by $|A| = (|a_{ij}|)$. The spectral radius of a square matrix $A$ is denoted by $\rho(A)$. A matrix $A = (a_{ij}) \in \mathbb{R}^{n\times n}$ is called a $Z$-matrix if $a_{ij}\leq 0$ for $i\neq j$ and $i, j = 1, 2, ..., n$, an $M$-matrix if $A$ is a $Z$-matrix with $A^{-1}\geq 0$, and a $P$-matrix if all of its principal minors are positive (^[25]). A matrix $A = (a_{ij}) \in \mathbb{R}^{n\times n}$ is called an $H$-matrix if its comparison matrix $\langle A\rangle$ is an $M$-matrix, and an $H_{+}$-matrix if it is an $H$-matrix with all positive diagonal elements, where the comparison matrix $\langle A\rangle = (\langle a_{ij}\rangle)$ is defined by $\langle a_{ii}\rangle = |a_{ii}|$ and $\langle a_{ij}\rangle = -|a_{ij}|$ with $i\neq j$ for $i, j = 1, 2, ..., n$ (^[2]). For a given matrix $A\in \mathbb{R}^{n\times n}$, the splitting $A = F - G$ is called an $M$-splitting if $F$ is an $M-$matrix and $G\geq 0$ (^[26]). For a given vector $x\in \mathbb{R}^{n}$, the symbols $x_{+}$ and $x_{-}$ denote the vectors $x_{+} = \max\{0, x\}$ and $x_{-} = \max\{0, -x\}$, respectively. For $x$, $x_{+}$ and $x_{-}$, we have

$\bf{Lemma\; 2.1}$ ^[27,28] Let $A\in \mathbb{R}^{n\times n}$ be an $M$-matrix and $A = F-G$ be an $M$-splitting. Then

For the equation

where $\Omega$ is a given positive diagonal matrix, we have the following conclusion.

$\bf{Lemma\; 2.2}$ The solution of the linear complementarity problem (1.1) and the solution of equation (2.1) have the following relation:

(ⅰ) If $x^{*}$ is a solution of (1.1) and $x^{**} = Ax^{*}+q$, then $x^{***} = x^{*}-\Omega x^{**}$ is a solution of (2.1).

(ⅱ) If $x^{*}$ is a solution of (2.1), then $x^{*}_{+}$ is a solution of (1.1).

Proof. (ⅰ) Suppose $x^{*}$ is a solution of (1.1), then we have

Therefore, denoting $Ax^{*}+q$ by $x^{**}$, for positive diagonal matrix $\Omega$, we have

It follows that

Thus

That is, $x^{***}$ is a solution of (2.1).

(ⅱ) Suppose $x^{*}$ is a solution of (2.1), i.e.,

Notice that $x^{*} = x^{*}_{+}-x^{*}_{-}$, we have $x^{*}_{-} = \Omega(Ax^{*}_{+}+q)$, and it follows that $Ax^{*}_{+}+q = \Omega^{-1} x^{*}_{-}\geq 0$. Moreover,

Therefore, $x^{*}_{+}$ is a solution of (1.1).

From Lemma 2.2 (ⅱ), we know that the solution of (1.1) can be obtained by solving Eq (2.1).

3.   General fixed-point method

In this section, we first prove that the solution of equation (2.1) is unique when the system matrix is a $P$-matrix, then propose the general fixed-point (GFP) method for solving (1.1) based on (2.1), and discuss the convergence conditions.

$\bf{Theorem\; 3.1}$ If $A$ is a $P$-matrix, then for any positive diagonal matrix $\Omega$, (2.1) has a unique solution.

Proof. Since $A$ is a $P$-matrix, the linear complementarity problem (1.1) has a unique solution for any $q\in \mathbb{R}^{n}$ (^[29]). Therefore, from Lemma 2.2, we know that Eq (2.1) has solution(s). Suppose $y^{*}, \; z^{*}$ are solutions of (2.1), then

By Lemma 2.2, $y_{+}^*$ and $z_{+}^*$ are are solutions of (1.1). Since (1.1) has a unique solution, $y^{*}_{+} = z^{*}_{+}$. It follows that $y^{*} = z^{*}$.

Based on (2.1) and Lemma 2.2 (ⅱ), we get an iterative method for solving (1.1):

We state the algorithm for (3.1) as follows.

$\bf{Theorem\; 3.2}$ Assume $A$ is a $P$-matrix. Let $\left\{x^{(k)}\right\}_{k = 1}^{+\infty}$ be the sequence generated by (3.1) and $x^{*}$ be the solution of (2.1). If

then $\left\{x^{(k)}\right\}_{k = 1}^{+\infty}$ converges to $x^{*}$ for any initial vector $x^{(0)}\in \mathbb{R}^n$.

Proof. Since $x^*$ is the unique solution of (2.1), we have

Combining with Eq (3.1), we get

It follows that

Therefore, if $\rho(|I-\Omega A|) < 1$, then $\left\{x^{(k)}\right\}_{k = 1}^{+\infty}$ converges to $x^{*}$ for any initial vector $x^{(0)}\in \mathbb{R}^n$.

Since $\rho(A)\leq\|A\|$ for any matrix norm $\|\cdot\|$, we can get the following corollary easily from Theorem 3.2.

$\bf{Corollary\; 3.1}$ Suppose $A$ is a $P$-matrix. Let $\left\{x^{(k)}\right\}_{k = 1}^{+\infty}$ be the sequence generated by (3.1) and $x^{*}$ be the solution of (2.1). If

then $\left\{x^{(k)}\right\}_{k = 1}^{+\infty}$ converges to $x^{*}$ for any initial vector $x^{(0)}\in \mathbb{R}^n$, where $\|\cdot\|$ is any matrix norm.

In the following, we consider two special cases: $A$ is an $H_{+}$ matrix with $\Omega = \omega D^{-1}$, where $D = \mathrm{diag}(A)$, and $A$ is a symmetric positive definite matrix with $\Omega = \omega I$, where $I$ is the identity matrix.

$\bf{Theorem\; 3.3}$ Suppose $A$ is an $H_{+}$-matrix, $D = \mathrm{diag}(A)$ and $B = D-A$. Let $\Omega = \omega D^{-1}$ with $\omega > 0$ and $\left\{x^{(k)}\right\}_{k = 1}^{+\infty}$ be the sequence generated by (3.1) and $x^{*}$ be the solution of (2.1). If

then $\left\{x^{(k)}\right\}_{k = 1}^{+\infty}$ converges to $x^{*}$ for any initial vector $x^{(0)}\in \mathbb{R}^n$. Moreover, $\omega = 1$ is the optimal choice.

Proof. Since $A$ is an $H_{+}$-matrix, $D = \mathrm{diag}(A)$, and $A = D-B$, we have that $\rho(D^{-1}|B|) < 1$; see ^[21]. For $\Omega = \omega D^{-1}$ with $\omega > 0$, we have

It follows that

It can be easily seen from (3.4) that $\rho(|I-\Omega A|) < 1$ for $\omega \in (0, 1]$ and for $\omega > 1$, $\rho(|I-\Omega A|) < 1$ if and only if $\omega < \frac{2}{1+\rho(D^{-1}|B|)}$. Therefore, if $0 < \omega < \frac{2}{1+\rho(D^{-1}|B|)}$, $\left\{x^{(k)}\right\}_{k = 1}^{+\infty}$ converges to $x^{*}$ for any initial vector $x^{(0)}\in \mathbb{R}^n$. Moreover, it can be seen from (3.4) that when $\omega = 1$, $\rho(|I-\Omega A|) = \rho(D^{-1}|B|)$ is minimal. That is, $\omega = 1$ is the optimal choice.

$\bf{Theorem\; 3.4}$ Suppose $A$ is a symmetric positive definite matrix. Set $\Omega = \omega I$ with $\omega > 0$ and denote the smallest and the largest eigenvalues of $A$ by $\lambda_{\min}$ and $\lambda_{\max}$, respectively. Let $\left\{x^{(k)}\right\}_{k = 1}^{+\infty}$ be the sequence generated by (3.1) and $x^{*}$ be the solution of (2.1). If

then $\left\{x^{(k)}\right\}_{k = 1}^{+\infty}$ converges to $x^{*}$ for any initial vector $x^{(0)}\in \mathbb{R}^n$.

Proof. From the proof of Theorem 3.2, similarly, we can have

Then

where $\parallel\cdot\parallel_2$ is the spectral norm of matrix. So, we have a convergence condition of (3.1), that is

Since

we can solve

and

to obtain the convergence conditions of (3.1), that is $0 < \omega\leq\frac{2}{\lambda_{\min}+\lambda_{\max}}$ and $\frac{2}{\lambda_{\min}+\lambda_{\max}}\leq\omega < \frac{2}{\lambda_{\max}}$, which can be combined into

Thus, the conclusion is proved.

Let $A$ be split as $A = D-L-U$, where $D$, $-L$ and $-U$ are the diagonal, strictly lower and strictly upper triangular parts of $A$, respectively. We can derive another iterative method for solving (1.1) based on (2.1) by using the idea of Gauss-Seidel:

We state the algorithm for (3.5) as follows.

We call iterative methods (3.1) and (3.5) the general fixed-point method. We note that in ^[4], Shi, Yang, and Huang introduced a fixed-point method, which is a special case of (3.1) with $\Omega = \alpha I$ ($\alpha > 0$). In the following, we analyze the convergence of iterative method (3.5). In particular, we consider the convergence domain of $\Omega$ when $A$ is an $H_{+}$-matrix.

$\bf{Theorem\; 3.5}$ Suppose $A$ is a $P$-matrix, then the sequence $\left\{x^{(k)}\right\}_{k = 1}^{+\infty}$ generated by (3.5) converges to the unique solution $x^{*}$ of (2.1) for any initial vector $x^{(0)}\in \mathbb{R}^n$ if

Proof. Since $A$ is a $P$-matrix, equation (2.1) has a unique solution for any positive diagonal matrix $\Omega$. Based on (2.1) and $A = D-L-U$, we get

From the above formula and equation (3.5), we have

Thus

It follows that

Since $I-|\Omega L|$ is an $M$-matrix, i.e., $(I-|\Omega L|)^{-1}\geq 0$, we have

Therefore, if $\rho((I-|\Omega L|)^{-1}|I-\Omega(D-U)|) < 1$, then $\left\{x^{(k)}\right\}_{k = 1}^{+\infty}$ converges to $x^{*}$ for any initial vector $x^{(0)}\in \mathbb{R}^n$.

We now consider the convergence domain of $\Omega$ for iterative method (3.5) when $A$ is an $H_{+}$-matrix.

$\bf{Theorem\; 3.6}$ Suppose $A$ is an $H_{+}$-matrix and either of the following conditions holds:

(1) $0 < \Omega\leq D^{-1}$;

(2) $\Omega > D^{-1}$ and $2\Omega^{-1}-D-|B|$ is an $M$-matrix, where $B = L+U$.

Then the sequence $\left\{x^{(k)}\right\}_{k = 1}^{+\infty}$ generated by (3.5) converges to the unique solution of (2.1) for any initial vector $x^{(0)}\in \mathbb{R}^n$.

Proof. Since any $H_{+}$-matrix is a $P$-matrix (^[2]), Eq (2.1) has a unique solution. Consider the splitting

(1) When $0 < \Omega\leq D^{-1}$, $(I-|\Omega L|)-|I-\Omega(D-U)|$ is an $M$-splitting of the $M$-matrix $\Omega \langle A\rangle$, therefore, it follows from Lemma 2.1 that $\rho((I-|\Omega L|)^{-1}\cdot|I-\Omega(D-U)|) < 1$.

(2) When $\Omega > D^{-1}$, if $2\Omega^{-1}-D-(|L|+|U|)$ is an $M$-matrix, then the splitting $(I-|\Omega L|)-|I-\Omega(D-U)|$ is an $M$-splitting of the $M$-matrix $\Omega(2\Omega^{-1}-D-|B|)$, therefore $\rho((I-|\Omega L|)^{-1}\cdot|I-\Omega(D-U)|) < 1$.

Collecting (1) and (2), this theorem is established based on Theorem 3.5.

$\bf{Corollary\; 3.2}$ Suppose $A$ is an $H_{+}$-matrix and $\Omega = \omega D^{-1}$ with $\omega > 0$, then the sequence $\left\{x^{(k)}\right\}_{k = 1}^{+\infty}$ generated by (3.5) converges to the unique solution $x^{*}$ of (2.1) for any initial vector $x^{(0)}\in \mathbb{R}^n$ if

Proof. For $0 < \omega \leq 1$, we have $\Omega \leq D^{-1}$, that is, $\Omega$ satisfies the first condition of Theorem 3.6. When $\Omega = \omega D^{-1}$, the second condition of Theorem 3.6 can be represented as

That is,

Therefore, if $0 < \omega < \frac{2}{1+\rho(D^{-1}|B|)}$ holds, $\left\{x^{(k)}\right\}_{k = 1}^{+\infty}$ converges to $x^{*}$.

To end this section, we consider the LCP$(A, q)$ arising in American option pricing problem, where $A$ is a symmetric tridiagonal $M$-matrix:

with $\lambda, \theta > 0$; see for instance ^[4]. In the following, we derive the optimal value of the parameter $\omega$ for the case $\Omega = \omega D^{-1}$ under the sense of $1$-norm.

For iterative method (3.1), based on (3.4), we consider

Let $\mu = \lambda \theta$. By some computation, we get

It can be easily seen that

That is, $\omega = 1$ is an optimal parameter.

We now consider the value of $\omega$ for iterative method (3.5). Let $\alpha = \frac{\omega}{1+2 \mu}$, i.e., $\omega = \alpha (1+2\mu)$, then

Let $\nu = |1- \omega|$, and $a_i = (\alpha \mu)^{i}$, $i = 0, 1, 2\ldots, n$. Then $a_{i} a_j = a_{i+j}$ and it is easy to get

Based on (3.6) and (3.8), we define

Since $a_i > 0$ and $0 < \nu < 1$, it can be seen that

We consider a particular $\omega$, that is $\omega = 1$, then

Remark. Based on the above discussions and notice that

we expect that iterative method (3.5) converges faster than iterative method (3.1) when $\omega = 1$.

4.   Numerical experiments

In this section, we illustrate some examples. Since most of the projected methods involve parameters, and it is not easy to select a proper in practice. Meanwhile, the implicit equations related to the projected methods are different, it is difficult to compare the GFP method with other projected methods fairly. So we do not compare the GFP method with other projected methods except for the FP method. The modulus-based approaches for solving the LCP$(A, q)$ include many cases, we select the modulus-based SOR (MSOR) iteration method with the best cases as comparison. Besides, we use the GFP method to solve the LCP$(A, q)$ arising in American option pricing problem(^[4]). The number of iteration steps, the elapsed time and the norm of the residual vector are denoted by IT, CPU and RES, respectively. RES is defined as:

where $x^{(k)}_{+}$ is the $k$th approximate solution of (1.1). The iteration process stops if $\mathrm{RES}(x^{(k)}_{+}) < 10^{-5}$ or the number of iteration steps reaches 1000.

The system matrix $A$ in the first two examples is generated by

here, $\mu$, $\eta$ and $\zeta$ are given constants, $\widehat{A} = \mathrm{Tridiag}(-I, S, -I) \in \mathbb{R}^{n\times n}$ is a block-tridiagonal matrix,

are two tridiagonal matrices, and $C = \mathrm{diag}([1\; 2\; 1\; 2\; \cdots])$ is a diagonal matrix of order $n$, $m$ and $n$ satisfy $n = m^{2}$. For convenience, we set

then the LCP($A(\mu, \eta, \zeta), q$) have a unique solution when $A(\mu, \eta, \zeta)$ is a $P$-matrix. All computations are run by using Matlab version 2016 on a Dell Laptop (Intel(R) Core(TM) i5-7200U CPU @ 2.50GHz, 4.00GB RAM).

$\bf{Example\; 4.1}$ In this example, we compare the GFP method with the MSOR iteration method (^[2]), which is a very effective method in solving the LCP($A, q$). Both the GFP method and the MSOR iteration method involve a parameter matrix $\Omega$, which has many choices and it is not appropriate if we take the same $\Omega$. We set $\Omega = \omega D^{-1}$ in the GFP method and $\Omega = \omega D$ in the MSOR iteration method, respectively. For fairness, we compare the two methods with the best cases by performing many experiments. Set $x^{(0)} = \mathrm{zeros}(n, 1)$ and $n = 10000$, then we obtain the following Table 1.

In Table 1, Alg 1 and Alg 2 denote Algorithm 1 and Algorithm 2, respectively. From Table 1, we can find that iterative method (3.1) is better than the MSOR iteration method in the running time, and iterative method (3.5) is better than the MSOR iteration method in the iteration steps. Since Algorithm 2 calculates the numerical solution entry by entry, the matrix-vector fast calculation technique can not be used, which makes the method more time-consuming. For Algorithm 1, single-step running time is short, but the number of iteration steps is relatively large. Besides, the best parameters for the MSOR iteration method are selected by many experiments, however, the best parameter for iterative method (3.1) is near or equals to 1, which is a conclusion of Theorem 3.3.

$\bf{Example\; 4.2}$ In this example, we illustrate the convergence domains of $\omega$ and the convergence rates of iterative methods (3.1) and (3.5). We consider $A(0, 0, 1)$ and $A(0, 1, 1)$. The former is a symmetric $H_{+}$-matrix and the latter is a nonsymmetric $H_{+}$-matrix, both of which are $P$-matrices. We set $\Omega = \omega D^{-1} = \omega\mathrm{diag}(A(\mu, \eta, \zeta))^{-1}$. From Theorem 3.3 and Corollary 3.2, we know that $\omega\in\left(0, \frac{2}{1+\rho(D^{-1}|B|)}\right)$ is a sufficient convergence domain for both iterative methods (3.1) and (3.5). Here, we consider a larger domain and set $\omega$ to be

where $\delta = \frac{1}{2}\left(\frac{2}{1+\rho(D^{-1}|B|)}-1\right)$. The symbols '$\mathrm{floor}$' and ':' are an integer function and a command in Matlab software, respectively. Then the right boundary of interval $\left(0, \frac{2}{1+\rho(D^{-1}|B|)}\right)$ is the penultimate one in these points. We denote $\rho(|I-\Omega A|)$ in (3.2) and $\rho((I-\Omega|L|)^{-1}|I-\Omega(D-U)|)$ in (3.6) by $\rho_1$ and $\rho_2$, respectively. Set $n = 900$ and $x^{(0)} = \mathrm{zeros}(n, 1)$, then we obtain Table 2 and the corresponding Figure 1 as follows.

From Table 2 and Figure 1, we can find that Theorem 3.3 and Corollary 3.2 provide the convergent domains of $\omega$ for iterative methods (3.1) and (3.5), respectively, and the initial iteration vector is arbitrary when $\omega$ falls in these domains. When $\omega$ exceeds the convergence domain, the two iterative methods are convergent sometimes when $x^{(0)} = \mathrm{zeros}(n, 1)$. Meanwhile, we can find that iterative method (3.5) is usually faster than iterative method (3.1) in terms of IT when $\omega$ takes the same value. In addition, this example illustrates the optimal parameter conclusion in Theorem 3.3, i.e., $\omega = 1$ is a good parameter for iterative method (3.1).

$\bf{Example\; 4.3}$ In this example, we apply the GFP method to solve the LCP ($A, q$) in ^[4], where the matrix $A$ satisfies (3.9). We set $\Omega = \omega D^{-1}$ and $\lambda\theta = 0.5, 1, 1.5, 2$, respectively. Similarly, just as Example 4.2, we set $\omega = 1-\mathrm{floor}\left(\frac{1}{\delta}\right)*\delta:\delta:\frac{2}{1+\rho(D^{-1}|B|)}+\delta$ with $\delta = \frac{1}{2}\left(\frac{2}{1+\rho(D^{-1}|B|)}-1\right)$, $\rho_{1} = \rho(|I-\Omega A|)$ in (3.2) and $\rho_{2} = \rho((I-\Omega|L|)^{-1}|I-\Omega(D-U)|)$ in (3.6). Then $\omega = 1$ is the fourth from the bottom of these values of $\omega$. We set $n = 900$ and $x^{(0)} = \mathrm{randn}(n, 1)$ in our experiments, then we obtain Figure 2 as follows.

From Figure 2, we can find that both iterative method (3.1) and iterative method (3.5) have good performance when $\omega = 1$, i.e., the two iterative methods solve the LCP ($A, q$) in very few iteration steps. Meanwhile, this example also verifies that iterative method (3.5) is faster than iterative method (3.1) when $\omega = 1$, which is a conclusion given at the end of Section 3.

$\bf{Example\; 4.4}$ In this example, we compare the GFP method with the FP method (^[4]). The system matrix is generated by

where $C = \mathrm{diag}([1\; 2\; \cdots\; n])$ and $n = 900$. The initial iteration vector is $x^{(0)} = \mathrm{zeros}(n, 1)$. For the FP method, the parameter matrix is $\Omega = \omega I$, and for the GFP method, the parameter matrix is $\Omega = \omega D^{-1} = \omega (\mathrm{diag}(A))^{-1}$. We consider two cases in our experiments, that is $\eta = 0$ and $\eta = 1$. Thus, these two matrices are $H_{+}$-matrices. Since the convergence domains of $\omega$ are different, we can not use the same $\omega$ values for the two methods. For the GFP method, based on Theorem 3.3 and Corollary 3.2, we know that the convergence domain is $0 < \omega < \frac{2}{1+\rho(D^{-1}|B|)}$, and we set $\omega$ to be $\frac{1}{3(1+\rho(D^{-1}|B|))}:\frac{1}{3(1+\rho(D^{-1}|B|))}:\frac{2}{1+\rho(D^{-1}|B|)}.$ For the FP method, there are two different situations. For $\eta = 0$, since the system matrix is a symmetric positive definite matrix, based on Theorem 3.4, we know that the convergence domain is $0 < \omega < \frac{2}{\lambda_{\max}}$, we set $\omega$ to be

For $\eta = 1$, by performing many experiments, we set $\omega$ to be $0.0003:0.0002:0.0013$, which include the convergence parameter values. The numerical results are shown in Tables 3 and 4, respectively.

From Tables 3 and 4, we can find that since $\rho$ is very large, the FP method can not obtain an approximate solution even the number of iteration steps reaches 1000. On the contrary, for the GFP method, both (3.1) and (3.5) can obtain the approximate solution in fewer iteration steps. Therefore, it is obvious that the GFP method is better than the FP method.

5.   Conclusions

In this paper, based on an equivalent fixed-point equation with a parameter matrix $\Omega$, we present the general fixed-point (GFP) method for solving the LCP($A, q$), which is the generalization of the fixed-point (FP) method. For this method, we discussed two iterative forms: one is the basic form and the other is the converted form, which is associated with matrix splitting. Both iterative forms can keep the spare structure of $A$ in the iteration processes, thus the sparse structure of $A$ can be applied to improve the effectiveness of this method. For the GFP method, the convergence conditions are proved and some concrete convergence domains of $\Omega$ as well as the optimal cases are presented. The iteration form of the GFP method is simple and the convergence rate is affected by the spectral radius of the iteration matrix. The numerical experiments show that the GFP method is an effective and competitive iterative method.

Acknowledgements

The author thanks the anonymous referees for providing many useful comments and suggestions that made this paper more readable. This work was supported by Zhaoqing Education and Development Project (No. ZQJYY2020093), the Characteristic Innovation Project of Department of Education of Guangdong Province (No. 2020KTSCX159), the Innovative Research Team Project of Zhaoqing University, the Scientific Research Ability Enhancement Program for Excellent Young Teachers of Zhaoqing University and Zhaoqing University Research Project (No. 611-612279).

Conflict of interest

The author confirms that there has no conflict of interest.