A shrinkage-thresholding projection method for sparsest solutions of LCPs

Abstract

In this paper, we study the sparsest solutions of linear complementarity problems (LCPs), which study has many applications, such as bimatrix games and portfolio selections. Mathematically, the underlying model is NP-hard in general. By transforming the complementarity constraints into a fixed point equation with projection type, we propose an $l_1$ regularization projection minimization model for relaxation. Through developing a thresholding representation of solutions for a key subproblem of this regularization model, we design a shrinkage-thresholding projection (STP) algorithm to solve this model and also analyze convergence of STP algorithm. Numerical results demonstrate that the STP method can efficiently solve this regularized model and get a sparsest solution of LCP with high quality.

MSC:90C33, 90C26, 90C90.

1 Introduction

Given a matrix $M\in R^{n\times n}$ and a vector $q\in R^n$, the linear complementarity problem, denoted by $LCP(q,M)$, is to find a vector $x\in R^n$ such that

$$x\ge 0,Mx+q\ge 0,x^T(Mx+q)=0.$$

The set of solutions to this problem is denoted by $SOL(q,M)$. Throughout the paper, we always suppose $SOL(q,M)\ne \mathrm{\varnothing }$.

Many real-world phenomena in engineering, physics, mechanics, and economics are governed by linear complementarity problems. Extensive studies of LCP have been done, see the books [1–3] and the references therein. Numerical methods for solving LCPs, such as the Newton method, the interior point method, and the nonsmooth equation method, have been extensively investigated in the literature. However, it seems that there is no study of the sparsest solutions for LCPs. In fact, in real applications, it is very necessary to investigate the sparsest solution of the LCPs. For example this is so in bimatrix games [1] and portfolio selections [4]. For more details, see [5].

In this paper, we try to find the sparsest solution of the LCP, which has the smallest number of nonzero entries. To be specific, we seek a vector $x\in R^n$ by solving the $l_0$-norm minimization problem,

$$\begin{array}{r}min{\parallel x\parallel }_0\\ \text{s.t.}x\ge 0,Mx+q\ge 0,x^T(Mx+q)=0,\end{array}$$

(1)

where ${\parallel x\parallel }_0$ stands for the number of nonzero components of x. A solution of (1) is called the sparsest solution of the LCP.

The above minimization problem (1) is in fact a sparse optimization with equilibrium constraints. From the point of view of constraint conditions, it is not easy to get solutions due to the equilibrium constraints, as well as the discontinuous objective function.

To overcome the difficulty, we make use of the ${\mathbf{F}}_{min}(x)$ C-functions to construct the penalty of violating the equilibrium constraints. Recall that a function $\psi :R^2\to R$ is called a C-function, where C stands for complementarity, if for any pair $(a,b)\in R^2$,

$$\psi (a,b)=0\iff a\ge 0,b\ge 0,ab=0.$$

Define ${\mathbf{F}}_{\psi }(x):R^n\to R^n$ by

$${\mathbf{F}}_{\psi }(x)=\left(\begin{array}{c}\psi (x_1,F_1(x))\\ ⋮\\ \psi (x_n,F_n(x))\end{array}\right).$$

The C-function ${\mathbf{F}}_{min}$, the extension of the ‘min’ function, is defined as follows:

$${\mathbf{F}}_{min}(x)=x-{\mathrm{\Pi }}_{R_{+}^n}(x-F(x)):=x-{[x-F(x)]}_{+},$$

(2)

where $F(x)=Mx+q$, with the vector $q\in R^n$ and the matrix $M\in R^{n\times n}$, and ${\mathrm{\Pi }}_{R_{+}^n}$ is the ‘Euclidean projector’ onto the nonnegative orthant. It is clear that

$$x\text{ solves the }LCP(q,M)\iff {\mathbf{F}}_{min}(x)=0.$$

(3)

Combining (2) and (3), we can obtain the equivalence between $LCP(q,M)$ and the fixed point equation ${\mathbf{F}}_{min}(x)=0$, that is,

$$x\text{ solves the }LCP(q,M)\iff x={[x-F(x)]}_{+}.$$

(4)

In the view of the objective function, problem (1) is an $l_0$-norm minimization problem, which is combination and generally NP-hard. The complexity of this model is generally proportional with the number of variables. In order to overcome such a difficultly, many researchers have suggested to relax the $l_0$ norm and instead consider the $l_1$ norm; see [6–11]. Hence, in this paper we consider applying the $l_1$ norm to find the sparsest solution of LCPs, and we obtain the following minimization problem to approximate (1):

$$\begin{array}{r}\underset{x\in R^n}{min}{\parallel x\parallel }_1\\ \text{s.t.}x={[x-F(x)]}_{+},\end{array}$$

(5)

where $F(x)=Mx+q$, ${\parallel x\parallel }_1={\sum }_{i=1}^n|x_i|$.

By further employing the $l_1$ regularization term for seeking sparsity, we introduce a new variable, $z\in R^n$, in order to simplify the objective function and carry on an alternative iteration. We apply the following $l_1$ regularization projection minimization problem to approximate (1):

$$\begin{array}{r}\underset{x,z\in R^n}{min}{f}_{\lambda }(x,z):={\parallel x-z\parallel }^2+\lambda {\parallel x\parallel }_1\\ \text{s.t.}z={[x-F(x)]}_{+},\end{array}$$

(6)

where $\lambda \in (0,\mathrm{\infty })$ is a given regularization parameter and $\parallel \cdot \parallel$ refers to the Euclidean norm. We call (6) the $l_1$ regularization projection minimization problem.

This paper is organized as follows. In Section 2, we approximate (1) by the $l_1$ regularization projection minimization problem (6), and show theoretically that (6) is a good approximation. In Section 3, we develop a shrinkage-thresholding representation theory for the subproblem of (6) and propose a shrinkage-thresholding projection (STP) algorithm for (6). The convergency of the STP algorithm is proved. Numerical results are demonstrated in Section 4 to show that (6) is promising to provide a sparsest solution of LCP.

2 The $l_1$ regularized approximation

In this section, we study the relation between the $l_1$ regularization projection model (6) and the original model (5), which indicates that the regularized model is a good approximation.

Theorem 2.1 For any fixed $\lambda >0$, the solution set of (6) is nonempty and bounded. Let $\{(x_{{\lambda }_k},z_{{\lambda }_k})\}$ be a solution sequence of (6), and $\{{\lambda }_k\}$ be any positive sequence converging to 0. If $SOL(q,M)\ne \mathrm{\varnothing }$, then $\{x_{{\lambda }_k}\}$ has at least one accumulation point, and any accumulation point $x^{\ast }$ of $\{x_{{\lambda }_k}\}$ is a solution of (5).

Proof For any fixed $\lambda >0$, it is easy to show the coercivity of the objective function $f_{\lambda }(x,z)$ in (6), which refers to the property that

$$f_{\lambda }(x,z)\to +\mathrm{\infty }\text{as }\parallel x\parallel \to \mathrm{\infty }\text{ and }\parallel z\parallel \to \mathrm{\infty }.$$

(7)

Also notice for all $x\in R^n$ and $z\in R^n$, $f_{\lambda }(x,z)\ge 0$. This, together with (7), implies that the level set

$$\mathcal{L}=\{(x,z)|f_{\lambda }(x,z)\le f_{\lambda }(x_0,z_0)\text{ and }z={[x-F(x)]}_{+}\}$$

is nonempty and compact, where $x_0\in R^n$ and $z_0={[x_0-F(x_0)]}_{+}$ are given points. The solution set of (6) is nonempty and bounded since $f_{\lambda }(x,z)$ is continuous on ℒ.

Now we show the second part of this theorem. Let $\hat{x}\in SOL(q,M)$ and $\hat{z}={[\hat{x}-F(\hat{x})]}_{+}$. From (4), we have $\hat{x}=\hat{z}$. Note that $(x_{{\lambda }_k},z_{{\lambda }_k})$ is a solution of (6) with $\lambda ={\lambda }_k$, and $z_{{\lambda }_k}={[x_{{\lambda }_k}-F(x_{{\lambda }_k})]}_{+}$. We have

$$\begin{array}{rl}max\{{\parallel x_{{\lambda }_k}-z_{{\lambda }_k}\parallel }^2,{\lambda }_k{\parallel x_{{\lambda }_k}\parallel }_1\}& \le {\parallel x_{{\lambda }_k}-z_{{\lambda }_k}\parallel }^2+{\lambda }_k{\parallel x_{{\lambda }_k}\parallel }_1\\ \le {\parallel \hat{x}-\hat{z}\parallel }^2+{\lambda }_k{\parallel \hat{x}\parallel }_1\\ ={\lambda }_k{\parallel \hat{x}\parallel }_1.\end{array}$$

(8)

Employing (8), we easily find that for any ${\lambda }_k>0$,

$${\parallel x_{{\lambda }_k}\parallel }_1\le {\parallel \hat{x}\parallel }_1.$$

(9)

Hence the whole sequence $\{x_{{\lambda }_k}\}$ is bounded, and it has at least one accumulation point. One easily notices that the sequence $\{z_{{\lambda }_k}\}$ is the same case as $\{x_{{\lambda }_k}\}$, which refers to the inequality ${\parallel x_{{\lambda }_k}-z_{{\lambda }_k}\parallel }^2\le {\lambda }_k{\parallel \hat{x}\parallel }_1$. Let $x^{\ast }$, $z^{\ast }$ be arbitrary accumulation points of $\{x_{{\lambda }_k}\}$ and $\{z_{{\lambda }_k}\}$, respectively, and $z_{{\lambda }_k}={[x_{{\lambda }_k}-F(x_{{\lambda }_k})]}_{+}$. Then there exists a subsequence of $\{{\lambda }_k\}$, say $\{{\lambda }_{k_j}\}$, such that

$$\underset{k_j\to \mathrm{\infty }}{lim}{x}_{{\lambda }_{k_j}}=x^{\ast }\text{and}\underset{k_j\to \mathrm{\infty }}{lim}{z}_{{\lambda }_{k_j}}=z^{\ast }.$$

We obtain $z^{\ast }={[x^{\ast }-F(x^{\ast })]}_{+}$ by letting $k_j$ tend to ∞ in $z_{{\lambda }_{k_j}}={[x_{{\lambda }_{k_j}}-F(x_{{\lambda }_{k_j}})]}_{+}$. Taking ${\lambda }_{k_j}\to 0$ to both sides of the inequality

$${\parallel x_{{\lambda }_{k_j}}-z_{{\lambda }_{k_j}}\parallel }^2\le {\lambda }_{k_j}{\parallel \hat{x}\parallel }_1,$$

we get $x^{\ast }=z^{\ast }$, which implies $x^{\ast }={[x^{\ast }-F(x^{\ast })]}_{+}$, that is, $x^{\ast }\in SOL(q,M)$. From ${\parallel x_{{\lambda }_{k_j}}\parallel }_1\le {\parallel \hat{x}\parallel }_1$ with $k_j$ tending to ∞, we get ${\parallel x^{\ast }\parallel }_1\le {\parallel \hat{x}\parallel }_1$. Then by the arbitrariness of $\hat{x}\in SOL(q,M)$, we know $x^{\ast }$ is a solution of problem (5). This completes the proof. □

3 Solution representation, algorithm, and convergence

Fixing $z^k\in R_{+}^n$, we consider an unconstrained minimization subproblem:

$$\underset{x\in R^n}{min}{f}_{\lambda }(x,z^k):={\parallel x-z^k\parallel }^2+\lambda {\parallel x\parallel }_1.$$

(10)

For fixed λ, a minimizer $x^{\ast }$ for the convex function (10) must satisfy the corresponding optimality conditions

$$x^{\ast }=S_{\lambda }\left(z^k\right),$$

(11)

where the shrinkage-thresholding operator $S_{\lambda }$ is defined by

$${\left(S_{\lambda }\left(z^k\right)\right)}_i=\{\begin{array}{cc}{z}_i^k-\frac{\lambda }{2},\hfill & z_i^k\ge \frac{\lambda }{2},\hfill \\ 0,\hfill & 0\le z_i^k<\frac{\lambda }{2}.\hfill \end{array}$$

(12)

It demonstrates that a solution $x^{\ast }\in R^n$ of the subproblem (10) can be analytically expressed by (11).

By the solution representation, we construct the following shrinkage-thresholding projection (STP) algorithm to solve the $l_1$ regularization projection minimization problem (6).

STP Algorithm

Now we will show that the STP algorithm is well defined, that is, (13) is implementable. Before doing this, we need the following lemmas.

Lemma 3.1 [12]

Let $P_{\mathrm{\Omega }}$ be a metric projection operator onto a nonempty closed convex set $\mathrm{\Omega }\in R^n$. Given $x\in R^n$ and $d\in R^n$, define

$$H(\alpha )=P_{\mathrm{\Omega }}[x-\alpha d],\alpha \ge 0,$$

then $\parallel H(\alpha )-x\parallel$ is nondecreasing with respect to α.

Lemma 3.2 The step size ${\alpha }_{k+1}$ in (13) must exist.

Proof From Lemma 3.1, we see that $\parallel x^{k+1}-{[x^{k+1}-\alpha F(x^{k+1})]}_{+}\parallel$ is nondecreasing with respect to α. It is obvious that $\alpha ({\parallel x^{k+1}-x^k\parallel }^2+{\parallel x^k-z^k\parallel }^2)$ is strictly increasing with respect to α since ${\parallel x^{k+1}-x^k\parallel }^2+{\parallel x^k-z^k\parallel }^2>0$ before the iterations stop. It follows that the term

$$g(\alpha ):={\parallel x^{k+1}-{[x^{k+1}-\alpha F\left(x^{k+1}\right)]}_{+}\parallel }^2+\alpha ({\parallel x^{k+1}-x^k\parallel }^2+{\parallel x^k-z^k\parallel }^2)$$

(13)

is strictly increasing with respect to α before the iterations stop. Thus $g(\beta {\gamma }^m)$ is strictly decreasing with respect to the nonnegative integer m before the iterations stop. Note that $x^{k+1}=S_{{\lambda }_k}(z^k)\in R_{+}^n$ and ${[x^{k+1}]}_{+}=x^{k+1}$, then we have

$$g(0)=\parallel x^{k+1}-{\left[x^{k+1}\right]}_{+}\parallel =0.$$

One notes that $S_{\lambda }(z^k)<z^k$ for $\lambda >0$ and any $z^k>0$. It follows that ${\parallel x^{k+1}-z^k\parallel }^2>0$ before the iterations stop. If $g(\beta {\gamma }^0)<{\parallel x^{k+1}-z^k\parallel }^2$, then

$${\alpha }_{k+1}=\beta$$

is just what we seek; if $g(\beta {\gamma }^0)\ge {\parallel x^{k+1}-z^k\parallel }^2$, there is a positive integer m such that

$$g\left(\beta {\gamma }^{m-1}\right)\ge {\parallel x^{k+1}-z^k\parallel }^2\text{and}g\left(\beta {\gamma }^m\right)<{\parallel x^{k+1}-z^k\parallel }^2,$$

then

$${\alpha }_{k+1}=\beta {\gamma }^m$$

is just what we seek. All in all, the step size ${\alpha }_{k+1}$ in (13) must exist. □

We now begin to analyze the convergence of the proposed STP algorithm.

Theorem 3.1 Let $\{(x^k,z^k)\}$ be a sequence generated by STP algorithm, then

• (i) $\{f_{{\lambda }_k}(x^k,z^k)\}$ is monotonically decreasing and converges to a constant $C^{\ast }$;

• (ii) $\{x^k\}$ and $\{z^k\}$ are bounded and suppose ${inf}_k{\alpha }_k=\alpha >0$, then $\{x^k\}$ and $\{z^k\}$ are both asymptotically regular, i.e.,

  $$\underset{k\to \mathrm{\infty }}{lim}\parallel x^{k+1}-x^k\parallel =0\mathit{\text{and}}\underset{k\to \mathrm{\infty }}{lim}\parallel z^{k+1}-z^k\parallel =0.$$

Moreover, any accumulation point of $\{x^k\}$ is a solution of $LCP(q,M)$.

Proof (i) From (10) and (11), we have

$$x^{k+1}=\underset{x\in R^n}{arg\hspace{0.17em}min}{f}_{{\lambda }_k}(x,z^k).$$

It follows that

$$f_{{\lambda }_k}(x^{k+1},z^k)\le f_{{\lambda }_k}(x^k,z^k).$$

(14)

Since ${\lambda }_k$ is monotonically decreasing, together with the inequality (13), we get

$$\begin{array}{r}{f}_{{\lambda }_{k+1}}(x^{k+1},z^{k+1})-f_{{\lambda }_k}(x^{k+1},z^k)\\ ={\parallel x^{k+1}-z^{k+1}\parallel }^2-{\parallel x^{k+1}-z^k\parallel }^2+({\lambda }_{k+1}-{\lambda }_k){\parallel x^{k+1}\parallel }_1\\ \le {\parallel x^{k+1}-{[x^{k+1}-{\alpha }_{k+1}F\left(x^{k+1}\right)]}_{+}\parallel }^2-{\parallel x^{k+1}-z^k\parallel }^2\\ <-{\alpha }_{k+1}({\parallel x^{k+1}-x^k\parallel }^2+{\parallel x^k-z^k\parallel }^2).\end{array}$$

Hence,

$$f_{{\lambda }_{k+1}}(x^{k+1},z^{k+1})<f_{{\lambda }_k}(x^{k+1},z^k)-{\alpha }_{k+1}({\parallel x^{k+1}-x^k\parallel }^2+{\parallel x^k-z^k\parallel }^2).$$

(15)

Combining (16) with (17), we get

$$\begin{array}{rcl}{f}_{{\lambda }_{k+1}}(x^{k+1},z^{k+1})& <& f_{{\lambda }_k}(x^k,z^k)-{\alpha }_{k+1}({\parallel x^{k+1}-x^k\parallel }^2+{\parallel x^k-z^k\parallel }^2)\\ \le & f_{{\lambda }_k}(x^k,z^k),\end{array}$$

(16)

which shows that $\{f_{{\lambda }_k}(x^k,z^k)\}$ is monotonically decreasing. Since $\{f_{{\lambda }_k}(x^k,z^k)\}$ is bounded from below, $\{f_{{\lambda }_k}(x^k,z^k)\}$ converges to a constant $C^{\ast }$. This verifies (i) of Theorem 3.1.

(ii) From the fact that $(x^k,z^k,{\lambda }^k)\in \{(x,z,\lambda )\in R^n\times R_{+}^n\times R_{+}:f(x,z,\lambda )\le f(x^0,z^0,{\lambda }_0)\}$, which is bounded, we see that $\{x^k\}$ and $\{z^k\}$ are both bounded.

From (18) and ${inf}_k{\alpha }_k=\alpha >0$, we have

$$\begin{array}{rcl}{f}_{{\lambda }_k}(x^k,z^k)-f_{{\lambda }_{k+1}}(x^{k+1},z^{k+1})& >& {\alpha }_{k+1}({\parallel x^{k+1}-x^k\parallel }^2+{\parallel x^k-z^k\parallel }^2)\\ \ge & \alpha ({\parallel x^{k+1}-x^k\parallel }^2+{\parallel x^k-z^k\parallel }^2).\end{array}$$

This then implies

$$\begin{array}{rcl}\sum _{k=0}^N({\parallel x^{k+1}-x^k\parallel }^2+{\parallel x^k-z^k\parallel }^2)& <& \frac{1}{\alpha }\sum _{k=0}^N[f_{{\lambda }_k}(x^k,z^k)-f_{{\lambda }_{k+1}}(x^{k+1},z^{k+1})]\\ =& \frac{1}{\alpha }[f_{{\lambda }_0}(x^0,z^0)-f_{{\lambda }_{N+1}}(x^{N+1},z^{N+1})]\\ \le & \frac{1}{\alpha }{f}_{{\lambda }_0}(x^0,z^0).\end{array}$$

Thus ${\sum }_{k=0}^{\mathrm{\infty }}({\parallel x^{k+1}-x^k\parallel }^2+{\parallel x^k-z^k\parallel }^2)$ is convergent, which yields

$$\parallel x^{k+1}-x^k\parallel \to 0\text{and}\parallel x^k-z^k\parallel \to 0\text{as }k\to +\mathrm{\infty }.$$

(17)

The above two limitations and the inequality

$$\parallel z^{k+1}-z^k\parallel \le \parallel z^{k+1}-x^{k+1}\parallel +\parallel x^{k+1}-x^k\parallel +\parallel x^k-z^k\parallel$$

yield

$$\parallel z^{k+1}-z^k\parallel \to 0\text{as }k\to +\mathrm{\infty }.$$

Since $\{x^k\}$ is bounded, $\{x^k\}$ has at least one accumulation. Let $x^{\ast }$ be an accumulation point of $\{x^k\}$ and a subsequence $\{x^{k_j}\}$ converging to $x^{\ast }$. Since $\{{\alpha }_{k_j}\}\subset [\alpha ,\beta ]$, without loss of generality, we suppose ${\alpha }_{k_j}\to \overline{\alpha }\in [\alpha ,\beta ]$ as $k_j\to \mathrm{\infty }$. It follows that

$$z^{k_j}={[x^{k_j}-{\alpha }_{k_j}F\left(x^{k_j}\right)]}_{+}\to z^{\ast }:={[x^{\ast }-\overline{\alpha }F\left(x^{\ast }\right)]}_{+}\text{as }{k}_j\to \mathrm{\infty }.$$

(18)

Combining (20) with (19), we get $x^{\ast }=z^{\ast }$, which gives $x^{\ast }={[x^{\ast }-\overline{\alpha }F(x^{\ast })]}_{+}$ and this means $x^{\ast }\in SOL(q,M)$. The proof is thus complete. □

4 Numerical experiments

In this section, we present some numerical experiments to demonstrate the effectiveness of our STP algorithm. All the numerical experiments were performed on a DELL computer (2.99 GHz, 2 GB of RAM), using MATLAB 7.10.

In the STP algorithm, the maximum number of iterations is set to 500. We end the STP algorithm, if $\parallel x^k-z^k\parallel <1.0\text{E}-5$ or if it reaches the maximum number of iterations. We set ${\lambda }_0=10$, $K=10$, $\tau =1/7$, $x^0$, and $z^0$ to be zero vectors as the initial points.

Test 1: Z-matrix LCPs [5]

Let us consider $LCP(q,M)$ where

$$M=I_n-\frac{1}{n}e{e}^T=\left(\begin{array}{cccc}1-\frac{1}{n}& -\frac{1}{n}& \cdots & -\frac{1}{n}\\ -\frac{1}{n}& 1-\frac{1}{n}& \cdots & -\frac{1}{n}\\ ⋮& ⋮& \ddots & ⋮\\ -\frac{1}{n}& -\frac{1}{n}& \cdots & 1-\frac{1}{n}\end{array}\right)\text{and}q=\left(\begin{array}{c}\frac{1}{n}-1\\ \frac{1}{n}\\ ⋮\\ \frac{1}{n}\end{array}\right).$$

Here $I_n$ is the identity matrix of order n and $e={(1,1,\dots ,1)}^T\in R^n$.

Such a matrix M is widely used in statistics. It is clear that M is a Z-matrix as well as a positive semi-definite matrix. For any scalars $a\ge 0$, we know that the vectors $x=ae+e_1$ are solutions to $LCP(q,M)$, because

$$x\ge 0,Mx+q=M{e}_1+q=0,x^T(Mx+q)=0.$$

Among all the solutions, the vector $x^{\ast }=e_1={(1,0,\dots ,0)}^T$ is the unique sparsest solution.

We test the STP algorithm for different dimensions with $n=100,500,1,000,3,000,5,000,7,000$, respectively. In this set of experiments, we set $\beta =1$ and $\gamma =1$. The results are displayed in Table 1.

Table 1 Computational results on LCPs with Z -matrices

In Table 1, ‘$D_{xz}$’ denotes the Euclidean distance between $x^k$ and $z^k$, which is in fact the value of the merit function at $x^k$, ‘gap’ denotes the Euclidean distance between $x^k$ and the true sparsest solution $x^{\ast }$, ‘Spar’ denotes the number of the entries of $x^k$ such that $x_i^k>1.0\text{E}-5$, and ‘time’ denotes the computational time in seconds.

From Table 1, we can see that the STP algorithm is effective to find the sparse solutions of LCPs. The sparsity of our solution $x^k$ is the same as the sparsity of the true sparse solution $x^{\ast }$. Moreover, the $D_{xz}$ and the ‘gap’ decrease as the dimension of the matrix M increases, which indicates that the larger size of the problem, the more effective the algorithm is.

Test 2: Randomly created LCPs with positive semidefinite matrices

In this subsection, we test STP for randomly created LCPs with positive semidefinite matrices.

First, we state the way of constructing LCPs and their solutions. Let a matrix $Z\in R^{n\times r}$ ($r<n$) be generated with the standard normal distribution and let $M=Z{Z}^T$. Let the sparse vector $\overline{x}$ be generated with the standard normal distribution and its sparsity be as follows: the sparsity is set to be $n/20$ if $n\le 1,000$; the sparsity is set to be $n/100$ if $1,000<n\le 5,000$; the sparsity is set to be $n/200$ if $n>5,000$. After the matrix M and the sparse vector $\overline{x}$ have been generated, a vector $q\in R^n$ can be constructed such that $\overline{x}$ is a solution of the $LCP(q,M)$. Then $\overline{x}$ can be regard as a sparse solution of the $LCP(q,M)$. Let M and q be input to our STP algorithm, then STP will output a solution $x^k$. We must emphasize that the sparsity of $x^k$ may be smaller than that of $\overline{x}$ since $\overline{x}$ maybe not the sparsest solution of the $LCP(q,M)$. In this case, $x^k$ is a sparser than $\overline{x}$.

In this set of experiments, ‘iter’ denotes the number of iterations for outputting $x^k$, ‘Spar-i’ denotes the number of the entry of $\overline{x}$ satisfying ${\overline{x}}_i^k>1.0\text{E}-5$, and ‘Spar-o’ denotes the number of the entry of $x^k$ satisfying $x_i^k>1.0\text{E}-5$. We set $\beta =0.9$ and $\gamma =0.5$.

We test the STP algorithm for different dimensions with $n=100,300,500,1,000,5,000,8,000$, respectively. For each case, we randomly run 10 times and compute the average values of ‘iter’, ‘Spar-i’, ‘Spar-o’ and $D_{xz}$, respectively. The results are displayed in Table 2.

Table 2 Results on randomly created LCPs with positive semidefinite matrices

From Table 2, we can see that the STP algorithm works very fast. Even for $n=8,000$, it only takes 29.7 seconds to yield a very sparse solution to the LCP. The values of $D_{xz}$ are all less than $1.0\text{E}-5$, which indicates the output points are solutions of $LCP(q,M)$. Moreover, the output solution $x^k$ is sparser than $\overline{x}$. When the dimension of M increases, the accuracy does not decrease but increases and the time cost by STP increases slowly. These phenomena show that STP is very robust. We can draw the conclusion that STP is very efficient for finding the sparsest solution of LCPs.

Remark The continuation method of the regularized parameter λ plays an important role in STP for find sparsest solutions of high quality. Moveover, a large amount of numerical experiments indicate that STP is very robust whenever $\lambda =1,5,10,20$.

5 Conclusions

In this paper, we concentrate on finding the sparsest solutions of LCPs. We propose an $l_1$ regularized projection minimization model. Then we develop a thresholding representation theory for the subproblem of $l_1$ regularized projection minimization problem, and design a shrinkage-thresholding projection (STP) algorithm to solve the regularized model. The convergence of the STP algorithm is proved. Preliminary numerical results indicate that the $l_1$ regularized model as well as the STP method are promising to find sparsest solutions of LCPs.