- Email: [email protected]
Conjugate gradient method for the linear complementarity problem with S -matrix
Mathematical and Computer Modelling 48 (2008) 918–928 www.elsevier.com/locate/mcm
Conjugate gradient method for the linear complementarity problem with S-matrixI Dong-Hui Li a,b,∗ , Yi-Yong Nie a , Jin-Ping Zeng b,c , Qing-Na Li b a Shenyang Institute of Automation, Chinese Academy of Sciences, China b College of Mathematics and Econometrics, Hunan University, Changsha, 410082, China c College of Software, Dongguan University of Technology, China
Received 16 November 2005; received in revised form 30 September 2007; accepted 23 October 2007
Abstract In this paper, we present a conjugate gradient method for solving the linear complementarity problem that involves an S-matrix. At each step, we solve a lower-dimensional system of linear equations by conjugate gradient method. The method terminates at the exact solution of the problem after a finite number of iterations. Moreover, the computational complexity of the proposed method is no more than the computational complexity of a conjugate gradient method for solving a system of linear equations. Preliminary numerical experiments show that the method is efficient. c 2007 Elsevier Ltd. All rights reserved.
Keywords: Linear complementarity problem; S-matrix; Conjugate gradient method
1. Introduction The conjugate gradient method is very efficient for solving systems of linear equations and nonlinear programming. It is particularly useful for solving large scale problems. When solving a system of linear equations or minimizing a strict convex quadratic function, the method terminates at the solution of the problem in finite steps. The purpose of this paper is to develop a conjugate gradient method for solving the linear complementarity problem. Let M = (m i j ) ∈ R n×n and q ∈ R n be given. The linear complementarity problem LCP(M, q) is to find an x ∈ R n such that M x − q ≥ 0,
x ≥ 0,
x T (M x − q) = 0.
Let N = {1, 2, . . . , n}. For x ∈ R n , define the index sets I (x) = {i ∈ N | xi = 0},
J (x) = N \ I (x).
I Supported by the NSF of China via grant 10771057 and the Doctoral Fund of Ministry of Education of China [2003] 172. ∗ Corresponding author at: Shenyang Institute of Automation, Chinese Academy of Sciences, China.
E-mail addresses: [email protected], [email protected] (D.-H. Li), [email protected] (J.-P. Zeng), [email protected] (Q.-N. Li). c 2007 Elsevier Ltd. All rights reserved. 0895-7177/$ - see front matter doi:10.1016/j.mcm.2007.10.017
(1.1)
D.-H. Li et al. / Mathematical and Computer Modelling 48 (2008) 918–928
919
Suppose that x ∗ is a solution of LCP(M, q). Let m T1 , m T2 , . . . , and m Tn denote the rows of M. It is clear that J (x ∗ ) = {i | xi∗ > 0} and that x ∗ is a solution of the following lower-dimensional system of linear equations: T m i x − qi = 0, i ∈ J (x ∗ ), xi = 0, i ∈ I (x ∗ ). On the other hand, if a solution x¯ of the system of linear equations T m i x − qi = 0, i ∈ J (x), xi = 0, i ∈ I (x)
(1.2)
satisfies x¯i ≥ 0,
∀i ∈ J (x) ¯
and
m iT x¯ − qi ≥ 0,
∀i ∈ I (x), ¯
then x¯ is a solution of LCP(M, q). In this case, we can solve the system of linear equations (1.2) to get a solution of the LCP(M, q). In general, solving a system of linear equations is much easier than solving an LCP. Therefore, it is interesting to find a system of linear equations whose solution is also a solution of the LCP. By using the so-called identification function [6], the LCP(M, q) can be locally converted into a system of linear equations. However, it is generally not easy to convert the LCP(M, q) into a system of linear equations like (1.2) because the index sets I (x ∗ ) and J (x ∗ ) depend on the solution x ∗ . In this paper, we consider the linear complementarity problem LCP(M, q) where M is an S-matrix, i.e. a symmetrical M-matrix. Here we use the definition of M-matrix given in [13]. In other word, Matrix M ∈ R n×n is called an M-matrix if it satisfies (i) M is nonsingular. (ii) M −1 is nonnegative. (iii) M is a Z -matrix, i.e. the elements of M satisfies m i j ≤ 0,
∀i 6= j, i, j = 1, 2, . . . , n.
If M is an S-matrix, it is also positive definite (see in [13]). Then it follows that LCP(M, q) is equivalent to the bound constrained strongly convex QP: 1 min x T M x − q T x. x≥0 2
(1.3)
Many large-scale scientific and engineering computation problems such as the finite element or finite difference discretization of free boundary problems are LCPs with S-matrices [3]. For such an LCP, we construct a finite sequence of systems of linear equations approaching the system of linear equations (1.2). The last system in this sequence is the system of linear equations (1.2). Specially, we first give a way to identify the index set J (x ∗ ). We do it by constructing a finite sequence of index sets {Jk }tk=0 satisfying J0 ⊂ J1 ⊂ · · · ⊂ Jt = J (x ∗ ). As a result, the LCP(M, q) is converted into the system of linear equations (1.2). For each k, 1 ≤ k ≤ t, the set Jk depends on the solution of the lower-dimensional system of linear equations ( m iT x − qi = 0, i ∈ Jk−1 (P(k − 1)) 4 xi = 0, i ∈ Ik−1 = N \ Jk−1 . We then propose a conjugate gradient method for solving LCP(M, q). The method is implemented by solving a sequence of lower-dimensional systems of linear equations P(1), P(2), . . . , P(t). For each k, problem P(k) is solved by a conjugate gradient method. Note that Jk−1 ⊂ Jk , k = 1, 2, . . . , t. The coefficient matrix M Jk−1 Jk−1 of problem (P(k − 1)) is a principal submatrix of the coefficient matrix M Jk Jk of problem P(k). When we solve problem P(k) by conjugate gradient method, we reuse the information obtained for solving problem (P(k − 1)) such that the total computation cost for solving the sequence of problem P(1), P(2), . . . , P(k) is no more than the computation cost of a conjugate gradient method for solving problem P(k). Consequently, the total computation cost of the proposed
920
D.-H. Li et al. / Mathematical and Computer Modelling 48 (2008) 918–928
conjugate gradient method for solving LCP(M, q) is no more than the computation cost of a conjugate gradient method for solving the lower-dimensional system of linear equations (1.2). The above conjugate gradient method is a direct method if no round-off error exists. Direct methods for solving the linear complementarity problem have received much attention in the last decades. We refer to [4] and [11] for comprehensive reviews in this field. For examples, the well-known Lemke’s method, the variable dimension method due to Van der Heyden, and the principal pivoting method due to Zoutendijk and Bard, etc. Usually, these direct methods are exponential methods (see, e.g. in [1,8,14]). There are also direct methods with polynomial computational complexity. Early in 1976, Mangasarian converted the linear complementarity problem into a linear program in some important special cases including those when M or its inverse is a Z -matrix [10]. As a result, it is solvable in polynomial time. Todd [16] obtained a polynomial probabilistic average computational complexity of the lexicographic Lemke method. The leading term of the complexity is 14 n 4 . In [17,18], direct methods for solving the linear complementarity with S-matrix and M-matrix were developed. The leading term of the computational complexity of these methods is n 4 . Early in [2], Chandrasekaran proposed a Lemke method. It was proved by Saigal [15] that when applied to solve a linear complementarity problem with Z -matrix, the leading term of the computational complexity of Chandrasekaran’s method is n 3 . Cryer [5] proposed an important modification of Saigal’s algorithm. When used to solve a linear complementarity problem with tridiagonal matrix, the complexity of the method is only O(n). Recently, by using a Gaussian eliminating technique, the authors in [9] developed a direct method for solving the linear complementarity problem with Z -matrix. The leading term of the computational complexity of the method is 13 n 3 . The aim of this paper is to develop a conjugate gradient method for solving the linear complementarity problem with S-matrix. It is particularly useful for solving large-scale problems. Moreover, as we shall show (see Section 3) that the computational complexity of the proposed method is no more than a conjugate gradient method for solving an n-dimensional system of linear equations. In the next section, we provide a way to estimate the index set J (x ∗ ). In Section 3, we propose the method and estimate its computational complexity. We conclude this section by introducing some notations used in the paper. We use E to denote the unit matrix whose dimension is clear from the context. Let x = (x1 , x2 , . . . , xn )T ∈ R n . For index set I ⊂ N , we use x I to denote the subvector of x whose elements are xi , i ∈ I . For matrix A = (ai j ) ∈ R n×n , and index sets I, J ⊂ N , we use A I J to denote the submatrix of A whose elements are ai j , i ∈ I , j ∈ J . For a finite set I ⊆ N , we use |I | to denote the number of elements in I . 2. Identification of positive variables In this section, we identify the positive variables of LCP(M, q) that forms the index set J (x ∗ ), where x ∗ is the solution of LCP(M, q). Throughout this section, we assume that M is an M-matrix. It is well known that LCP(M, q) has a unique solution if M is an M-matrix (see, e.g. in [4,11]). Define the index sets J0 = {i | qi > 0},
I0 = N \ I0 = {i | qi ≤ 0}.
(2.1)
It is clear that if J0 = ∅, then the solution of LCP(M, q) is x ∗ = 0. Suppose that J0 6= ∅. As we mentioned in Section 1, the following lemma is obvious but very important: Lemma 2.1. The unique solution of LCP(M, q) coincides with the unique solution of the following lowerdimensional system of linear equations: T m i x − qi = 0, i ∈ J (x ∗ ), (2.2) xi = 0, i ∈ I (x ∗ ). On the other hand, if for some index sets I, J ⊆ N satisfying I ∩ J = ∅ and I ∪ J = N , the solution of the lower-dimensional system of linear equations T m i x − qi = 0, i ∈ J, (2.3) xi = 0, i ∈ I
D.-H. Li et al. / Mathematical and Computer Modelling 48 (2008) 918–928
921
satisfies x J ≥ 0,
and
M I J x J − q I ≥ 0,
(2.4)
then x is the solution of LCP(M, q). The above lemma establishes the equivalence between LCP(M, q) and the system of linear equations. If we can identify the index J (x ∗ ), then LCP(M, q) can be converted into a lower-dimensional system of linear equations. In what follows, we give a way to identify the index set J (x ∗ ). We fulfil this purpose by constructing a monotone sequence of subsets of J (x ∗ ). The idea was originated in [18]. First, we introduce a lemma from [9]. For completeness we include the proof in the Appendix. Lemma 2.2. Let the index sets α, β satisfy α ⊆ J (x ∗ ) = {i ∈ N | xi∗ > 0}, α ∩ β = ∅ and α ∪ β = N . Denote E αα (Mαα )−1 Mαβ (Mαα )−1 qα ¯ M= , q¯ = , (2.5) 0 Mββ − Mβα (Mαα )−1 Mαβ qβ − Mβα (Mαα )−1 qα where E αα ∈ R |α|×|α| is an identity matrix. Then the unique solution of LCP(M, q) coincides with the unique solution ¯ q). of LCP( M, ¯ The following lemma shows that J0 is a subset of J (x ∗ ). Lemma 2.3. We have J0 ⊆ J (x ∗ ). Proof. Since x ∗ is the solution LCP(M, q), we have for each i ∈ N , X m ii xi∗ ≥ − m i j x ∗j + qi ≥ qi . j6=i
For any i ∈ J0 , qi > 0. Therefore, we must have xi∗ > 0. This means i ∈ J (x ∗ ).
Let x 0 be the solution of P(0). If condition (2.4) is satisfied with I = I0 and J = J0 , then it follows form Lemma 2.1 that J0 = J (x ∗ ) and x 0 is the solution of LCP(M, q). In this case, we have identified the index set J (x ∗ ). Suppose that (2.4) does not hold with x = x 0 and I = I0 , J = J0 . We are going to find a larger index set J1 such that J0 ⊂ J1 ⊆ J (x ∗ ). The next lemma gives a way to find such index J1 . Lemma 2.4. Let x 0 be the unique solution of P(0) and (2.4) be not satisfied with I = I0 and J = J0 . Define the index set J10 = {i ∈ I0 | m iT x 0 − qi < 0}. Then we have J10 6= ∅ and J10 ⊂ J (x ∗ ). Proof. The fact J10 = 6 ∅ is obvious. We turn to verifying J10 ⊂ J (x ∗ ). Let E J0 J0 (M J0 J0 )−1 M J0 I0 (M J0 J0 )−1 q J0 ¯ M= , q¯ = . 0 M I0 I0 − M I0 J0 (M J0 J0 )−1 M J0 I0 q I0 − M I0 J0 (M J0 J0 )−1 q J0
(2.6)
¯ q) It follows from Lemma 2.2 that LCP(M, q) and LCP( M, ¯ are equivalent. Note that i ∈ J10 if and only if i ∈ I0 and q¯i > 0. Therefore, we claim from Lemma 2.3 that i ∈ J (x ∗ ). It is obvious that J10 ∩ J0 = ∅. Letting J1 = J0 ∪ J10 , we get a new index set J1 which contains J0 as a proper subset. Moreover, it follows from Lemma 2.4 that J1 ⊆ J (x ∗ ). Note that matrix M¯ defined by (2.6) is also an M¯ q) matrix and J1 consists of the index i for which q¯i > 0. Index set J1 in problem LCP( M, ¯ plays the same role as the index J0 in problem LCP(M, q). In particular, we have the following corollary:
922
D.-H. Li et al. / Mathematical and Computer Modelling 48 (2008) 918–928
Corollary 2.5. The solution of the lower-dimensional system of linear equations M J1 J1 x J1 − q J1 = 0 is positive. Repeating the above process, we get a sequence of index sets J0 ⊂ J1 ⊂ · · · ⊂ Jt ⊆ J (x ∗ ). Since J (x ∗ ) is finite, after finite steps, we get for some t ≤ n that Jt = J (x ∗ ). The subscript t is the integer for which the solution x t of problem P(t) satisfies (2.4) with J = Jt and I = It ≡ N \ Jt . The above process of generating {Jk }tk=0 essentially provides a method for solving LCP(M, q). It is the direct method developed by Zhou [18] for LCP(M, q) with S-matrix M. In the next section, by means of the identification of J (x ∗ ), we propose a conjugate gradient method for solving LCP(M, q). The computation cost of the proposed method is no more than the cost of a conjugate gradient method for solving the lower-dimensional system of linear equations (1.2). 3. The conjugate gradient algorithm In this section, we propose a conjugate gradient method for solving LCP(M, q). Throughout this section, we assume that matrix M is an S-matrix, i.e. it is a symmetrical M-matrix. It is clear that any S-matrix must be positive definite [13]. Therefore, the linear complementarity problem LCP(M, q) is equivalent to the following strictly convex quadratic programming: ( 4 1 min f (x) = x T M x − q T x, (3.1) 2 s.t. x ≥ 0. The idea of the proposed conjugate gradient method is to solve a sequence of lower-dimensional systems of linear equations by the conjugate gradient method. To save computation cost, we use the conjugate directions in the previous iteration to generate the conjugate directions for the next iteration. The total computation cost of the proposed method is no more than the cost of the conjugate gradient method for solving the lower-dimensional system of linear equations (1.2), which is no more than the computation cost of a conjugate gradient method for solving an n-dimensional system of linear equations. We first simply recall the process of a conjugate gradient method for solving the system of linear equations: Ax − b = 0,
(3.2)
where A ∈ R n×n is symmetrical and positive definite, and b ∈ R n . This system of linear equations is equivalent to the following unconstrained optimization problem: min f (x) =
1 T x Ax − bT x. 2
(3.3)
Let d 0 , d 1 , . . . , d n−1 be A-conjugate directions, i.e. (d i )T Ad j = 0,
∀i 6= j, i, j = 0, 1, . . . , n − 1.
Starting from an arbitrary x 0 ∈ R n , the conjugate gradient method generates a sequence of iterates {x j } by letting x j+1 = x j + α j d j ,
j = 0, 1, . . . ,
where α j is determined by αj = −
∇ f (x j )T d j (Ax j − b)T d j =− , j T j (d ) Ad (d j )T Ad j
j = 0, 1, . . . .
If we neglect the round-off errors, then the above process terminates at the solution of (3.2) within n steps. In other words, there is an integer t ≤ n such that x t = A−1 b.
D.-H. Li et al. / Mathematical and Computer Modelling 48 (2008) 918–928
923
We now use this idea to construct a conjugate gradient method for solving LCP(M, q). Let J0 ⊆ J (x ∗ ) be defined by (2.1). For k ≥ 0, let Jk ⊂ J (x ∗ ) be generated by the process described in Section 2 and i k be the number of elements in Jk . Let d J0k , d J1k , . . . , d iJkk−1 be M Jk Jk -conjugate directions in R ik . Starting from x 0Jk , we solve the system of linear equations M Jk Jk x Jk − q Jk = 0
(3.4)
by conjugate gradient method and get x 1Jk , x 2Jk , . . . , x iJkk . That is, j+1
x Jk
j
j
= x Jk + α j d Jk ,
j = 0, 1, . . . , i k − 1,
where j
j
αj = −
(M Jk Jk x Jk − q Jk )T d Jk j
j
(d Jk )T M Jk Jk d Jk
,
j = 0, 1, . . . , i k − 1.
Since M Jk Jk is an S-matrix, it is positive definite. By the property of the conjugate gradient method, in the case where no round-off errors exist, there is an index tk ≤ i k such that x tJkk is the solution of the system of linear equations (3.4). Moreover, it follows from Corollary 2.5 that x tJkk = (M Jk Jk )−1 q Jk > 0. Letting Ik = N \ Jk and extending x tJkk to x tk = x tJkk ∪Ik by setting x Itkk = 0 Ik , we get an estimate x tk of x ∗ . To save computation cost, we establish a connection between the k-th system of linear equations and the (k + 1)-th system. Note that Jk+1 = Jk ∪ Jk0 , where Jk0 satisfies Jk0 ⊂ J (x ∗ ) and Jk0 ∩ Jk = ∅. So, M Jk+1 Jk+1 can be partitioned as the following four blocks: M Jk Jk0 M Jk Jk . M Jk0 Jk M Jk0 Jk0 For i = 0, 1, . . . , tk − 1, extend d iJk to d iJk+1 by setting d iJk0 = 0 Jk0 ,
i = 0, 1, . . . , tk − 1.
Since d iJk , i = 0, . . . , tk − 1 are M Jk Jk -conjugate directions, it is not difficult to show that d iJk+1 , i = 0, . . . , tk − 1 are M Jk+1 Jk+1 -conjugate directions. Therefore, while solving the system of linear equations M Jk+1 Jk+1 x Jk+1 − q Jk+1 = 0
(3.5)
by conjugate gradient method, we have already had M Jk+1 Jk+1 -conjugate directions d iJk+1 , i = 0, . . . , tk − 1. Extend x tJkk to x tJkk+1 by setting x tJkk0 = 0 Jk0 ,
i = 0, 1, . . . , tk − 1.
It is easy to verify that x iJk+1 , i = 0, 1, 2, . . . , tk satisfy: j+1
j
j
x Jk+1 = x Jk+1 + α j d Jk+1 ,
j = 0, 1, . . . , tk − 1,
924
D.-H. Li et al. / Mathematical and Computer Modelling 48 (2008) 918–928
where j
j
αj = −
(M Jk Jk x Jk − b Jk )T d Jk j
j
(d Jk )T M Jk Jk d Jk
=−
(M Jk+1 Jk+1 x Jk+1 − q Jk+1 )T d Jk+1 j
j
(d Jk+1 )T M Jk+1 Jk+1 d Jk+1
,
j = 0, 1, . . . , i k − 1.
tk −1 In other words, the sequence {x iJk+1 }i=0 can be regarded as the iterates generated by the conjugate gradient method for solving the system of linear equations (3.5) with conjugate directions d iJk+1 , i = 0, 1, . . . , tk − 1. We summarize the above process as the algorithm below.
Algorithm 1 (Conjugate Gradient Method for LCP). 4
4
Step 0 Let J0 = {i | qi > 0} and I0 = N \ J0 . If J0 is empty, we let I0 = N = {1, . . . , n}. Given initial point x 0J0 and M J0 J0 -conjugate directions d iJ0 , i = 0, 1, . . . , |J0 | − 1. Let k := 0 and t−1 = 0. t
Step 1 Starting from x Jk−1 , solve the system of linear equations k M Jk Jk x Jk − q Jk = 0
(3.6)
by conjugate gradient method and get the solution x tJkk . Step 2 Set Ik = N \ Jk . If M Ik Jk x tJkk − q Ik ≥ 0, then we get the solution x ∗Jk
=
x tJkk ,
x∗
(3.7) of LCP(M, q) by setting
= 0 Ik .
x I∗k
Otherwise, go to Step 3. Step 3 Let ( ) X tk Jk0 = i ∈ Ik m i j x j − qi < 0 , j∈J
and
Jk+1 := Jk ∪ Jk0 .
k
Extend x tJkk and d iJk , i = 0, 1, . . . , tk − 1, to x tJkk+1 and d iJk+1 , i = 0, 1, . . . , tk − 1, respectively, by setting x tJkk0 = 0 Jk0 ,
d iJk0 = 0 Jk0 ,
i = 0, 1, . . . , tk − 1.
Step 4 Construct M Jk+1 Jk+1 -conjugate directions d iJk+1 , i = tk , tk + 1, . . . , |Jk+1 | − 1 such that directions d iJk+1 , i = 0, 1, . . . , |Jk+1 | − 1 are M Jk+1 Jk+1 -conjugate directions. Let k := k + 1. Go to Step 1. Remark. (i) Algorithm 1 can be viewed in some sense as an active set method for solving bound constrained strongly convex quadratic programming (1.3). The active set methods are very efficient for solving (1.3) (see e.g. [7] and references therein). However, Algorithm 1 has some advantages. As pointed out in the previous part of this section, to solve the system of linear equations (3.6), we only need to partly use the conjugate gradient method. That is, in Step 1 of Algorithm 1, while solving the lower-dimensional system of linear equations (3.6) by conjugate gradient method, we have already had points j
xk ,
j = 0, 1, . . . , tk−1 j
which were obtained from previous iterations, where xk denotes the i-th element of xk . We only need to proceed the t process, starting from x Jk−1 and doing line searches along directions k j
d Jk ,
j = tk−1 , tk−1 + 1, . . . , |Jk | − 1.
In other words, Algorithm 1 generates the sequence of iterates {x k } that takes the following form: ! ! t t x J00 x 0J0 x J00 x 0J0 = , . . . , = ,..., t 0 I0 x I00 x I00 0 I0
D.-H. Li et al. / Mathematical and Computer Modelling 48 (2008) 918–928
t
x Jj−1 j
+1
t
+1
j−1 t +1 = x J j j−1 0I j xI
t
! ,...,
x Jjj
!
t
j
x I jj
t
x Jjj = 0I j
925
! ,....
(ii) We did not give the way to generate the M Jk Jk -conjugate directions. It can be done by any existing conjugate gradient methods such as the Fletcher–Reeves method and the Polak–Ribi`ere method [12]. Algorithm 1 can be regarded as a conjugate gradient method for solving the lower-dimensional system of linear equations (1.2). By means of the well-known convergence of the conjugate gradient method, we immediately have the following convergence theorem (see e.g. [12, Theorem 5.5]). Theorem 3.1. Let M be an S-matrix. Then Algorithm 1 terminates at the solution of LCP(M, q) in at most |J (x ∗ )| steps. Moreover, we have the following convergence rate estimate: √ k √ λ1 − λt k+1 ∗ kx − x k M J (x ∗ )J (x ∗ ) ≤ 2 √ kx 0 − x ∗ k M J (x ∗ )J (x ∗ ) , √ λ1 + λt 4 p T where kxk M J (x ∗ )J (x ∗ ) = x M J (x ∗ )J (x ∗ ) x, and λ1 and λt are the largest and smallest eigenvalues of matrix M J (x ∗ )J (x ∗ ) . 4. Numerical experiments In this section, we present some numerical results to test the performance of the proposed method. The test problem comes from the finite difference discretization of the one side obstacle problem [3]: h−4u − f, v − ui ≥ 0, where K = {v ∈ as follows:
H01 (Ω )
M x − q ≥ 0,
∀v ∈ K ,
: v ≥ 0}, f = 4 sin(4x y), Ω = (0, 1) × (0, 1). By discretization, we obtain the test problem x T (M x − q) = 0,
x ≥ 0,
where h = r1 , n = r 2 , q = (4h 2 sin( 4ir 2j ))i j , i, j = 1, . . . , r , B −I M =
−I B .. .
−I .. . ..
.
.. ..
.
. −I
−I B
∈ R n×n ,
I is the identity matrix of r -dimension, and 4 −1 −1 B=
4 .. .
−1 .. . ..
.
.. ..
.
∈ R r ×r .
. −1 −1 4
We test the method on the problem with different dimensions and different initial points. The method was coded in Matlab 6.5 and run on a personal computer with 1.73 GHZ CPU processor. We stop the iteration if the condition k min{M x + q, x}k ≤ 10−7 is satisfied. At iteration k, we generate the M Jk Jk -conjugate direction by use of the standard FR method. We stop the inner iteration if the condition kM Jk Jk x Jk − q Jk k ≤ 10−5 is satisfied. Tables 1 and 2 list the results with different initial points x 0J0 = (1, . . . , 1)T and x 0J0 = (0, . . . , 0)T , respectively. The meaning of each column in these tables is as follows: “n” denotes the dimension of the problem; “|J0 |” stands for the number of elements in J0 ; “InIt” stands for the total number of inner iterations, namely, the total number of CG-iterations to
926
D.-H. Li et al. / Mathematical and Computer Modelling 48 (2008) 918–928
Table 1 Test results for Algorithm 1 with x0 = (1, . . . , 1)T n
InIt
OutIt
|J0 |
CPU
100 400 900 1600 2500 3600 4900 6400 8100 10 000 12 100 14 400 16 900 19 600 22 500
26 53 83 118 156 193 230 264 292 314 340 371 395 420 445
4 6 8 9 11 14 15 18 19 20 22 23 25 26 28
94 385 872 1551 2426 3497 4764 6222 7877 9730 11 775 14 014 16 451 19 084 21 906
0.063 0.172 0.343 0.703 1.313 2.234 3.688 5.515 7.938 11.187 14.594 19.766 24.938 30.672 37.375
Table 2 Test results for Algorithm 1 with x0 = (0, . . . , 0)T n
InIt
OutIt
|J0 |
CPU
100 400 900 1600 2500 3600 4900 6400 8100 10 000 12 100 14 400 16 900 19 600 22 500
20 40 63 93 125 157 187 213 237 256 285 305 323 338 364
4 6 8 9 11 14 15 17 18 20 22 23 25 26 28
94 385 872 1551 2426 3497 4764 6222 7877 9730 11 775 14 014 16 451 19 084 21 906
0.046 0.14 0.281 0.578 1.11 1.937 3.25 4.969 6.75 9.172 12.734 17 20.875 25.953 32.063
solve linear equations (3.6); “OutIt” means the total number of outer iterations and “CPU” denotes the computer time used (in second); We see from Tables 1 and 2 that starting from each initial point, the method terminates at the solution of the problem successfully. The proposed method can solve large-scale problems efficiently. 5. Discussion We have developed a conjugate gradient method for solving the linear complementarity problem with S-matrix. The method is an extension of the conjugate gradient method for solving systems of linear equations and minimization problems. The method reserves good properties of the conjugate gradient method for solving system of linear equations and minimization problems. An attractive property of the proposed method is that the computation cost of the method is no more than the cost of a conjugate gradient method for solving an n-dimensional system of linear equations. The proposed method depends on useful properties of S-matrix. These properties were used to identify the positive variables. It is interesting to develop a conjugate gradient method for solving the linear complementarity problem with general symmetrical and positive definite matrix. The difficulty to do this lies in the identification of the positive
D.-H. Li et al. / Mathematical and Computer Modelling 48 (2008) 918–928
927
variables. It is also interesting to develop preconditioned conjugate gradient method. This seems not trivial because Algorithm 1 implicitly implies the solution of a sequence of systems of linear equations. We leave it as a topic for further research. Acknowledgements The authors would like to thank the two anonymous referees for their very helpful comments. Appendix The Proof of Lemma 2.2. Since M is an M-matrix, its principal submatrices Mαα and Mββ are also M-matrices (see, e.g., in [13]), and then the Schur complement Mββ − Mβα (Mαα )−1 Mαβ is a Z -matrix. Noting that Iαα 0 Mαα Mαβ Iαα −(Mαα )−1 Mαβ Mβα Mββ −Mβα (Mαα )−1 Iββ 0 Iββ Mαα 0 = 0 Mββ − Mβα (Mαα )−1 Mαβ and Mββ − Mβα (Mαα )−1 Mαβ = Bββ , where
Bαα Bβα
Bαβ Bββ
=
Mαα Mβα
Mαβ Mββ
−1 ≥ 0,
we have that Mββ − Mβα (Mαα )−1 Mαβ is nonsingular and [Mββ − Mβα (Mαα )−1 Mαβ ]−1 ≥ 0. Hence, Mββ − Mβα (Mαα )−1 Mαβ is an M-matrix and then M¯ is an M-matrix by (2.5). Therefore, LCP(M, q) and ¯ q) LCP( M, ¯ exist unique solution respectively. So, it is sufficient to prove that the solution of LCP(M, q) is also the ¯ q). solution of LCP( M, ¯ ∗ Let x = x ≥ 0 be the solution of LCP(M, q). Since α ⊂ J (x ∗ ), by Lemma 2.1, we have m iT x − qi = 0,
i ∈ α,
i.e., Mαα xα + Mαβ xβ − qα = 0. Therefore, xα + (Mαα )−1 Mαβ xβ − (Mαα )−1 qα = 0.
(A.8)
And then (M x − q)β = Mβα xα + Mββ xβ − qβ = Mβα (Mαα )−1 (qα − Mαβ xβ ) + Mββ xβ − qβ = [Mββ − Mβα (Mαα )−1 Mαβ ]xβ − [qβ − Mβα (Mαα )−1 qα ]. Since x =
x∗
≥ 0 is the solution of LCP(M, q), by (A.8) and (A.9), we get
xα ≥ 0,
( M¯ x − q) ¯ α = xα + (Mαα )−1 Mαβ xβ − (Mαα )−1 qα = 0,
xβ ≥ 0,
( M¯ x − q) ¯ β = (M x − q)β ≥ 0,
and xβT ( M¯ x − q) ¯ β = 0,
¯ q). which implies x is the solution of LCP( M, ¯ The proof is then completed.
(A.9)
928
D.-H. Li et al. / Mathematical and Computer Modelling 48 (2008) 918–928
References [1] J.R. Birge, A. Gana, Computational complexity of Van der Heyden’s variable dimension algorithm and Dantzig–Cottle’s principal pivoting method for solving LCP’s, Mathematical Programming 26 (1983) 316–325. [2] R. Chandrasekaran, A special case of the complementarity pivot problem, Opsearch 7 (1970) 263–268. [3] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland Publishing Company, Amsterdam, 1978. [4] R.W. Cottle, J-S. Pang, R.E. Stone, The Linear Complementarity Problem, Academic Press, Inc., 1992. [5] C.W. Cryer, The eficient solution of linear complementarity problems for tradiagonal Minkowski matrices, ACM Transactions on Mathematical Software 9 (1983) 199–214. [6] F. Facchinei, A. Fischer, C. Kanzow, On the accurate identification of active constraints, SIAM Journal on Optimization 9 (1999) 14–32. [7] W.W. Hager, H. Zhang, A new active set algorithm for box constrained optimization, SIAM Journal on Optimization 17 (2007) 526–557. [8] H. Kremers, D. Talman, A new pivoting algorithm for the complementarity problem allowing for an arbitrary starting point, Mathematical Programming 63 (1994) 235–252. [9] D.H. Li, J.P. Zeng, Z.Z. Zhang, Gaussian pivoting method for solving linear complementarity problem, Applied Mathematics - A Journal of Chinese Universities 12 (1997) 419–426. [10] O.L. Mangasarian, Linear complementarity problems solvable by A single linear program, Mathematical Programming 10 (1976) 263–270. [11] K.G. Murty, Linear Complementarity, Linear and Nonlinear Programming, Heldermann, Berlin, 1988. [12] J. Nocedal, S. Wright, Numerical Optimization, Springer-Verlag, 2000. [13] J.M. Ortega, W.C. Rheinboldt, Iterative Solution of Nonlinear Equations in Several Variables, Academic Press, New York, 1970. [14] J. Rohn, A short proof of finiteness of Murty’s principal pivoting algorithm, Mathematical Programming 46 (1990) 255–256. [15] R. Saigal, Lemke’s algorithm and a special class of linear complementarity problems, Opsearch 8 (1971) 201–208. [16] M.J. Todd, Polynomial expected behavior of a pivoting algorithm for linear complementarity and linear programming problems, Mathematical Programming 35 (1986) 173–192. [17] L. Zhang, X.Y. Hu, On the direct method for linear complementarity problem, Mathematik Numerica Sinica 1 (1994) 59–64. [18] S.Z. Zhou, A direct method for the linear complementarity problem, Journal of Computational Mathematics 2 (1990) 178–182.
Conjugate gradient method for the linear complementarity problem with S-matrixI Dong-Hui Li a,b,∗ , Yi-Yong Nie a , Jin-Ping Zeng b,c , Qing-Na Li b a Shenyang Institute of Automation, Chinese Academy of Sciences, China b College of Mathematics and Econometrics, Hunan University, Changsha, 410082, China c College of Software, Dongguan University of Technology, China
Received 16 November 2005; received in revised form 30 September 2007; accepted 23 October 2007
Abstract In this paper, we present a conjugate gradient method for solving the linear complementarity problem that involves an S-matrix. At each step, we solve a lower-dimensional system of linear equations by conjugate gradient method. The method terminates at the exact solution of the problem after a finite number of iterations. Moreover, the computational complexity of the proposed method is no more than the computational complexity of a conjugate gradient method for solving a system of linear equations. Preliminary numerical experiments show that the method is efficient. c 2007 Elsevier Ltd. All rights reserved.
Keywords: Linear complementarity problem; S-matrix; Conjugate gradient method
1. Introduction The conjugate gradient method is very efficient for solving systems of linear equations and nonlinear programming. It is particularly useful for solving large scale problems. When solving a system of linear equations or minimizing a strict convex quadratic function, the method terminates at the solution of the problem in finite steps. The purpose of this paper is to develop a conjugate gradient method for solving the linear complementarity problem. Let M = (m i j ) ∈ R n×n and q ∈ R n be given. The linear complementarity problem LCP(M, q) is to find an x ∈ R n such that M x − q ≥ 0,
x ≥ 0,
x T (M x − q) = 0.
Let N = {1, 2, . . . , n}. For x ∈ R n , define the index sets I (x) = {i ∈ N | xi = 0},
J (x) = N \ I (x).
I Supported by the NSF of China via grant 10771057 and the Doctoral Fund of Ministry of Education of China [2003] 172. ∗ Corresponding author at: Shenyang Institute of Automation, Chinese Academy of Sciences, China.
E-mail addresses: [email protected], [email protected] (D.-H. Li), [email protected] (J.-P. Zeng), [email protected] (Q.-N. Li). c 2007 Elsevier Ltd. All rights reserved. 0895-7177/$ - see front matter doi:10.1016/j.mcm.2007.10.017
(1.1)
D.-H. Li et al. / Mathematical and Computer Modelling 48 (2008) 918–928
919
Suppose that x ∗ is a solution of LCP(M, q). Let m T1 , m T2 , . . . , and m Tn denote the rows of M. It is clear that J (x ∗ ) = {i | xi∗ > 0} and that x ∗ is a solution of the following lower-dimensional system of linear equations: T m i x − qi = 0, i ∈ J (x ∗ ), xi = 0, i ∈ I (x ∗ ). On the other hand, if a solution x¯ of the system of linear equations T m i x − qi = 0, i ∈ J (x), xi = 0, i ∈ I (x)
(1.2)
satisfies x¯i ≥ 0,
∀i ∈ J (x) ¯
and
m iT x¯ − qi ≥ 0,
∀i ∈ I (x), ¯
then x¯ is a solution of LCP(M, q). In this case, we can solve the system of linear equations (1.2) to get a solution of the LCP(M, q). In general, solving a system of linear equations is much easier than solving an LCP. Therefore, it is interesting to find a system of linear equations whose solution is also a solution of the LCP. By using the so-called identification function [6], the LCP(M, q) can be locally converted into a system of linear equations. However, it is generally not easy to convert the LCP(M, q) into a system of linear equations like (1.2) because the index sets I (x ∗ ) and J (x ∗ ) depend on the solution x ∗ . In this paper, we consider the linear complementarity problem LCP(M, q) where M is an S-matrix, i.e. a symmetrical M-matrix. Here we use the definition of M-matrix given in [13]. In other word, Matrix M ∈ R n×n is called an M-matrix if it satisfies (i) M is nonsingular. (ii) M −1 is nonnegative. (iii) M is a Z -matrix, i.e. the elements of M satisfies m i j ≤ 0,
∀i 6= j, i, j = 1, 2, . . . , n.
If M is an S-matrix, it is also positive definite (see in [13]). Then it follows that LCP(M, q) is equivalent to the bound constrained strongly convex QP: 1 min x T M x − q T x. x≥0 2
(1.3)
Many large-scale scientific and engineering computation problems such as the finite element or finite difference discretization of free boundary problems are LCPs with S-matrices [3]. For such an LCP, we construct a finite sequence of systems of linear equations approaching the system of linear equations (1.2). The last system in this sequence is the system of linear equations (1.2). Specially, we first give a way to identify the index set J (x ∗ ). We do it by constructing a finite sequence of index sets {Jk }tk=0 satisfying J0 ⊂ J1 ⊂ · · · ⊂ Jt = J (x ∗ ). As a result, the LCP(M, q) is converted into the system of linear equations (1.2). For each k, 1 ≤ k ≤ t, the set Jk depends on the solution of the lower-dimensional system of linear equations ( m iT x − qi = 0, i ∈ Jk−1 (P(k − 1)) 4 xi = 0, i ∈ Ik−1 = N \ Jk−1 . We then propose a conjugate gradient method for solving LCP(M, q). The method is implemented by solving a sequence of lower-dimensional systems of linear equations P(1), P(2), . . . , P(t). For each k, problem P(k) is solved by a conjugate gradient method. Note that Jk−1 ⊂ Jk , k = 1, 2, . . . , t. The coefficient matrix M Jk−1 Jk−1 of problem (P(k − 1)) is a principal submatrix of the coefficient matrix M Jk Jk of problem P(k). When we solve problem P(k) by conjugate gradient method, we reuse the information obtained for solving problem (P(k − 1)) such that the total computation cost for solving the sequence of problem P(1), P(2), . . . , P(k) is no more than the computation cost of a conjugate gradient method for solving problem P(k). Consequently, the total computation cost of the proposed
920
D.-H. Li et al. / Mathematical and Computer Modelling 48 (2008) 918–928
conjugate gradient method for solving LCP(M, q) is no more than the computation cost of a conjugate gradient method for solving the lower-dimensional system of linear equations (1.2). The above conjugate gradient method is a direct method if no round-off error exists. Direct methods for solving the linear complementarity problem have received much attention in the last decades. We refer to [4] and [11] for comprehensive reviews in this field. For examples, the well-known Lemke’s method, the variable dimension method due to Van der Heyden, and the principal pivoting method due to Zoutendijk and Bard, etc. Usually, these direct methods are exponential methods (see, e.g. in [1,8,14]). There are also direct methods with polynomial computational complexity. Early in 1976, Mangasarian converted the linear complementarity problem into a linear program in some important special cases including those when M or its inverse is a Z -matrix [10]. As a result, it is solvable in polynomial time. Todd [16] obtained a polynomial probabilistic average computational complexity of the lexicographic Lemke method. The leading term of the complexity is 14 n 4 . In [17,18], direct methods for solving the linear complementarity with S-matrix and M-matrix were developed. The leading term of the computational complexity of these methods is n 4 . Early in [2], Chandrasekaran proposed a Lemke method. It was proved by Saigal [15] that when applied to solve a linear complementarity problem with Z -matrix, the leading term of the computational complexity of Chandrasekaran’s method is n 3 . Cryer [5] proposed an important modification of Saigal’s algorithm. When used to solve a linear complementarity problem with tridiagonal matrix, the complexity of the method is only O(n). Recently, by using a Gaussian eliminating technique, the authors in [9] developed a direct method for solving the linear complementarity problem with Z -matrix. The leading term of the computational complexity of the method is 13 n 3 . The aim of this paper is to develop a conjugate gradient method for solving the linear complementarity problem with S-matrix. It is particularly useful for solving large-scale problems. Moreover, as we shall show (see Section 3) that the computational complexity of the proposed method is no more than a conjugate gradient method for solving an n-dimensional system of linear equations. In the next section, we provide a way to estimate the index set J (x ∗ ). In Section 3, we propose the method and estimate its computational complexity. We conclude this section by introducing some notations used in the paper. We use E to denote the unit matrix whose dimension is clear from the context. Let x = (x1 , x2 , . . . , xn )T ∈ R n . For index set I ⊂ N , we use x I to denote the subvector of x whose elements are xi , i ∈ I . For matrix A = (ai j ) ∈ R n×n , and index sets I, J ⊂ N , we use A I J to denote the submatrix of A whose elements are ai j , i ∈ I , j ∈ J . For a finite set I ⊆ N , we use |I | to denote the number of elements in I . 2. Identification of positive variables In this section, we identify the positive variables of LCP(M, q) that forms the index set J (x ∗ ), where x ∗ is the solution of LCP(M, q). Throughout this section, we assume that M is an M-matrix. It is well known that LCP(M, q) has a unique solution if M is an M-matrix (see, e.g. in [4,11]). Define the index sets J0 = {i | qi > 0},
I0 = N \ I0 = {i | qi ≤ 0}.
(2.1)
It is clear that if J0 = ∅, then the solution of LCP(M, q) is x ∗ = 0. Suppose that J0 6= ∅. As we mentioned in Section 1, the following lemma is obvious but very important: Lemma 2.1. The unique solution of LCP(M, q) coincides with the unique solution of the following lowerdimensional system of linear equations: T m i x − qi = 0, i ∈ J (x ∗ ), (2.2) xi = 0, i ∈ I (x ∗ ). On the other hand, if for some index sets I, J ⊆ N satisfying I ∩ J = ∅ and I ∪ J = N , the solution of the lower-dimensional system of linear equations T m i x − qi = 0, i ∈ J, (2.3) xi = 0, i ∈ I
D.-H. Li et al. / Mathematical and Computer Modelling 48 (2008) 918–928
921
satisfies x J ≥ 0,
and
M I J x J − q I ≥ 0,
(2.4)
then x is the solution of LCP(M, q). The above lemma establishes the equivalence between LCP(M, q) and the system of linear equations. If we can identify the index J (x ∗ ), then LCP(M, q) can be converted into a lower-dimensional system of linear equations. In what follows, we give a way to identify the index set J (x ∗ ). We fulfil this purpose by constructing a monotone sequence of subsets of J (x ∗ ). The idea was originated in [18]. First, we introduce a lemma from [9]. For completeness we include the proof in the Appendix. Lemma 2.2. Let the index sets α, β satisfy α ⊆ J (x ∗ ) = {i ∈ N | xi∗ > 0}, α ∩ β = ∅ and α ∪ β = N . Denote E αα (Mαα )−1 Mαβ (Mαα )−1 qα ¯ M= , q¯ = , (2.5) 0 Mββ − Mβα (Mαα )−1 Mαβ qβ − Mβα (Mαα )−1 qα where E αα ∈ R |α|×|α| is an identity matrix. Then the unique solution of LCP(M, q) coincides with the unique solution ¯ q). of LCP( M, ¯ The following lemma shows that J0 is a subset of J (x ∗ ). Lemma 2.3. We have J0 ⊆ J (x ∗ ). Proof. Since x ∗ is the solution LCP(M, q), we have for each i ∈ N , X m ii xi∗ ≥ − m i j x ∗j + qi ≥ qi . j6=i
For any i ∈ J0 , qi > 0. Therefore, we must have xi∗ > 0. This means i ∈ J (x ∗ ).
Let x 0 be the solution of P(0). If condition (2.4) is satisfied with I = I0 and J = J0 , then it follows form Lemma 2.1 that J0 = J (x ∗ ) and x 0 is the solution of LCP(M, q). In this case, we have identified the index set J (x ∗ ). Suppose that (2.4) does not hold with x = x 0 and I = I0 , J = J0 . We are going to find a larger index set J1 such that J0 ⊂ J1 ⊆ J (x ∗ ). The next lemma gives a way to find such index J1 . Lemma 2.4. Let x 0 be the unique solution of P(0) and (2.4) be not satisfied with I = I0 and J = J0 . Define the index set J10 = {i ∈ I0 | m iT x 0 − qi < 0}. Then we have J10 6= ∅ and J10 ⊂ J (x ∗ ). Proof. The fact J10 = 6 ∅ is obvious. We turn to verifying J10 ⊂ J (x ∗ ). Let E J0 J0 (M J0 J0 )−1 M J0 I0 (M J0 J0 )−1 q J0 ¯ M= , q¯ = . 0 M I0 I0 − M I0 J0 (M J0 J0 )−1 M J0 I0 q I0 − M I0 J0 (M J0 J0 )−1 q J0
(2.6)
¯ q) It follows from Lemma 2.2 that LCP(M, q) and LCP( M, ¯ are equivalent. Note that i ∈ J10 if and only if i ∈ I0 and q¯i > 0. Therefore, we claim from Lemma 2.3 that i ∈ J (x ∗ ). It is obvious that J10 ∩ J0 = ∅. Letting J1 = J0 ∪ J10 , we get a new index set J1 which contains J0 as a proper subset. Moreover, it follows from Lemma 2.4 that J1 ⊆ J (x ∗ ). Note that matrix M¯ defined by (2.6) is also an M¯ q) matrix and J1 consists of the index i for which q¯i > 0. Index set J1 in problem LCP( M, ¯ plays the same role as the index J0 in problem LCP(M, q). In particular, we have the following corollary:
922
D.-H. Li et al. / Mathematical and Computer Modelling 48 (2008) 918–928
Corollary 2.5. The solution of the lower-dimensional system of linear equations M J1 J1 x J1 − q J1 = 0 is positive. Repeating the above process, we get a sequence of index sets J0 ⊂ J1 ⊂ · · · ⊂ Jt ⊆ J (x ∗ ). Since J (x ∗ ) is finite, after finite steps, we get for some t ≤ n that Jt = J (x ∗ ). The subscript t is the integer for which the solution x t of problem P(t) satisfies (2.4) with J = Jt and I = It ≡ N \ Jt . The above process of generating {Jk }tk=0 essentially provides a method for solving LCP(M, q). It is the direct method developed by Zhou [18] for LCP(M, q) with S-matrix M. In the next section, by means of the identification of J (x ∗ ), we propose a conjugate gradient method for solving LCP(M, q). The computation cost of the proposed method is no more than the cost of a conjugate gradient method for solving the lower-dimensional system of linear equations (1.2). 3. The conjugate gradient algorithm In this section, we propose a conjugate gradient method for solving LCP(M, q). Throughout this section, we assume that matrix M is an S-matrix, i.e. it is a symmetrical M-matrix. It is clear that any S-matrix must be positive definite [13]. Therefore, the linear complementarity problem LCP(M, q) is equivalent to the following strictly convex quadratic programming: ( 4 1 min f (x) = x T M x − q T x, (3.1) 2 s.t. x ≥ 0. The idea of the proposed conjugate gradient method is to solve a sequence of lower-dimensional systems of linear equations by the conjugate gradient method. To save computation cost, we use the conjugate directions in the previous iteration to generate the conjugate directions for the next iteration. The total computation cost of the proposed method is no more than the cost of the conjugate gradient method for solving the lower-dimensional system of linear equations (1.2), which is no more than the computation cost of a conjugate gradient method for solving an n-dimensional system of linear equations. We first simply recall the process of a conjugate gradient method for solving the system of linear equations: Ax − b = 0,
(3.2)
where A ∈ R n×n is symmetrical and positive definite, and b ∈ R n . This system of linear equations is equivalent to the following unconstrained optimization problem: min f (x) =
1 T x Ax − bT x. 2
(3.3)
Let d 0 , d 1 , . . . , d n−1 be A-conjugate directions, i.e. (d i )T Ad j = 0,
∀i 6= j, i, j = 0, 1, . . . , n − 1.
Starting from an arbitrary x 0 ∈ R n , the conjugate gradient method generates a sequence of iterates {x j } by letting x j+1 = x j + α j d j ,
j = 0, 1, . . . ,
where α j is determined by αj = −
∇ f (x j )T d j (Ax j − b)T d j =− , j T j (d ) Ad (d j )T Ad j
j = 0, 1, . . . .
If we neglect the round-off errors, then the above process terminates at the solution of (3.2) within n steps. In other words, there is an integer t ≤ n such that x t = A−1 b.
D.-H. Li et al. / Mathematical and Computer Modelling 48 (2008) 918–928
923
We now use this idea to construct a conjugate gradient method for solving LCP(M, q). Let J0 ⊆ J (x ∗ ) be defined by (2.1). For k ≥ 0, let Jk ⊂ J (x ∗ ) be generated by the process described in Section 2 and i k be the number of elements in Jk . Let d J0k , d J1k , . . . , d iJkk−1 be M Jk Jk -conjugate directions in R ik . Starting from x 0Jk , we solve the system of linear equations M Jk Jk x Jk − q Jk = 0
(3.4)
by conjugate gradient method and get x 1Jk , x 2Jk , . . . , x iJkk . That is, j+1
x Jk
j
j
= x Jk + α j d Jk ,
j = 0, 1, . . . , i k − 1,
where j
j
αj = −
(M Jk Jk x Jk − q Jk )T d Jk j
j
(d Jk )T M Jk Jk d Jk
,
j = 0, 1, . . . , i k − 1.
Since M Jk Jk is an S-matrix, it is positive definite. By the property of the conjugate gradient method, in the case where no round-off errors exist, there is an index tk ≤ i k such that x tJkk is the solution of the system of linear equations (3.4). Moreover, it follows from Corollary 2.5 that x tJkk = (M Jk Jk )−1 q Jk > 0. Letting Ik = N \ Jk and extending x tJkk to x tk = x tJkk ∪Ik by setting x Itkk = 0 Ik , we get an estimate x tk of x ∗ . To save computation cost, we establish a connection between the k-th system of linear equations and the (k + 1)-th system. Note that Jk+1 = Jk ∪ Jk0 , where Jk0 satisfies Jk0 ⊂ J (x ∗ ) and Jk0 ∩ Jk = ∅. So, M Jk+1 Jk+1 can be partitioned as the following four blocks: M Jk Jk0 M Jk Jk . M Jk0 Jk M Jk0 Jk0 For i = 0, 1, . . . , tk − 1, extend d iJk to d iJk+1 by setting d iJk0 = 0 Jk0 ,
i = 0, 1, . . . , tk − 1.
Since d iJk , i = 0, . . . , tk − 1 are M Jk Jk -conjugate directions, it is not difficult to show that d iJk+1 , i = 0, . . . , tk − 1 are M Jk+1 Jk+1 -conjugate directions. Therefore, while solving the system of linear equations M Jk+1 Jk+1 x Jk+1 − q Jk+1 = 0
(3.5)
by conjugate gradient method, we have already had M Jk+1 Jk+1 -conjugate directions d iJk+1 , i = 0, . . . , tk − 1. Extend x tJkk to x tJkk+1 by setting x tJkk0 = 0 Jk0 ,
i = 0, 1, . . . , tk − 1.
It is easy to verify that x iJk+1 , i = 0, 1, 2, . . . , tk satisfy: j+1
j
j
x Jk+1 = x Jk+1 + α j d Jk+1 ,
j = 0, 1, . . . , tk − 1,
924
D.-H. Li et al. / Mathematical and Computer Modelling 48 (2008) 918–928
where j
j
αj = −
(M Jk Jk x Jk − b Jk )T d Jk j
j
(d Jk )T M Jk Jk d Jk
=−
(M Jk+1 Jk+1 x Jk+1 − q Jk+1 )T d Jk+1 j
j
(d Jk+1 )T M Jk+1 Jk+1 d Jk+1
,
j = 0, 1, . . . , i k − 1.
tk −1 In other words, the sequence {x iJk+1 }i=0 can be regarded as the iterates generated by the conjugate gradient method for solving the system of linear equations (3.5) with conjugate directions d iJk+1 , i = 0, 1, . . . , tk − 1. We summarize the above process as the algorithm below.
Algorithm 1 (Conjugate Gradient Method for LCP). 4
4
Step 0 Let J0 = {i | qi > 0} and I0 = N \ J0 . If J0 is empty, we let I0 = N = {1, . . . , n}. Given initial point x 0J0 and M J0 J0 -conjugate directions d iJ0 , i = 0, 1, . . . , |J0 | − 1. Let k := 0 and t−1 = 0. t
Step 1 Starting from x Jk−1 , solve the system of linear equations k M Jk Jk x Jk − q Jk = 0
(3.6)
by conjugate gradient method and get the solution x tJkk . Step 2 Set Ik = N \ Jk . If M Ik Jk x tJkk − q Ik ≥ 0, then we get the solution x ∗Jk
=
x tJkk ,
x∗
(3.7) of LCP(M, q) by setting
= 0 Ik .
x I∗k
Otherwise, go to Step 3. Step 3 Let ( ) X tk Jk0 = i ∈ Ik m i j x j − qi < 0 , j∈J
and
Jk+1 := Jk ∪ Jk0 .
k
Extend x tJkk and d iJk , i = 0, 1, . . . , tk − 1, to x tJkk+1 and d iJk+1 , i = 0, 1, . . . , tk − 1, respectively, by setting x tJkk0 = 0 Jk0 ,
d iJk0 = 0 Jk0 ,
i = 0, 1, . . . , tk − 1.
Step 4 Construct M Jk+1 Jk+1 -conjugate directions d iJk+1 , i = tk , tk + 1, . . . , |Jk+1 | − 1 such that directions d iJk+1 , i = 0, 1, . . . , |Jk+1 | − 1 are M Jk+1 Jk+1 -conjugate directions. Let k := k + 1. Go to Step 1. Remark. (i) Algorithm 1 can be viewed in some sense as an active set method for solving bound constrained strongly convex quadratic programming (1.3). The active set methods are very efficient for solving (1.3) (see e.g. [7] and references therein). However, Algorithm 1 has some advantages. As pointed out in the previous part of this section, to solve the system of linear equations (3.6), we only need to partly use the conjugate gradient method. That is, in Step 1 of Algorithm 1, while solving the lower-dimensional system of linear equations (3.6) by conjugate gradient method, we have already had points j
xk ,
j = 0, 1, . . . , tk−1 j
which were obtained from previous iterations, where xk denotes the i-th element of xk . We only need to proceed the t process, starting from x Jk−1 and doing line searches along directions k j
d Jk ,
j = tk−1 , tk−1 + 1, . . . , |Jk | − 1.
In other words, Algorithm 1 generates the sequence of iterates {x k } that takes the following form: ! ! t t x J00 x 0J0 x J00 x 0J0 = , . . . , = ,..., t 0 I0 x I00 x I00 0 I0
D.-H. Li et al. / Mathematical and Computer Modelling 48 (2008) 918–928
t
x Jj−1 j
+1
t
+1
j−1 t +1 = x J j j−1 0I j xI
t
! ,...,
x Jjj
!
t
j
x I jj
t
x Jjj = 0I j
925
! ,....
(ii) We did not give the way to generate the M Jk Jk -conjugate directions. It can be done by any existing conjugate gradient methods such as the Fletcher–Reeves method and the Polak–Ribi`ere method [12]. Algorithm 1 can be regarded as a conjugate gradient method for solving the lower-dimensional system of linear equations (1.2). By means of the well-known convergence of the conjugate gradient method, we immediately have the following convergence theorem (see e.g. [12, Theorem 5.5]). Theorem 3.1. Let M be an S-matrix. Then Algorithm 1 terminates at the solution of LCP(M, q) in at most |J (x ∗ )| steps. Moreover, we have the following convergence rate estimate: √ k √ λ1 − λt k+1 ∗ kx − x k M J (x ∗ )J (x ∗ ) ≤ 2 √ kx 0 − x ∗ k M J (x ∗ )J (x ∗ ) , √ λ1 + λt 4 p T where kxk M J (x ∗ )J (x ∗ ) = x M J (x ∗ )J (x ∗ ) x, and λ1 and λt are the largest and smallest eigenvalues of matrix M J (x ∗ )J (x ∗ ) . 4. Numerical experiments In this section, we present some numerical results to test the performance of the proposed method. The test problem comes from the finite difference discretization of the one side obstacle problem [3]: h−4u − f, v − ui ≥ 0, where K = {v ∈ as follows:
H01 (Ω )
M x − q ≥ 0,
∀v ∈ K ,
: v ≥ 0}, f = 4 sin(4x y), Ω = (0, 1) × (0, 1). By discretization, we obtain the test problem x T (M x − q) = 0,
x ≥ 0,
where h = r1 , n = r 2 , q = (4h 2 sin( 4ir 2j ))i j , i, j = 1, . . . , r , B −I M =
−I B .. .
−I .. . ..
.
.. ..
.
. −I
−I B
∈ R n×n ,
I is the identity matrix of r -dimension, and 4 −1 −1 B=
4 .. .
−1 .. . ..
.
.. ..
.
∈ R r ×r .
. −1 −1 4
We test the method on the problem with different dimensions and different initial points. The method was coded in Matlab 6.5 and run on a personal computer with 1.73 GHZ CPU processor. We stop the iteration if the condition k min{M x + q, x}k ≤ 10−7 is satisfied. At iteration k, we generate the M Jk Jk -conjugate direction by use of the standard FR method. We stop the inner iteration if the condition kM Jk Jk x Jk − q Jk k ≤ 10−5 is satisfied. Tables 1 and 2 list the results with different initial points x 0J0 = (1, . . . , 1)T and x 0J0 = (0, . . . , 0)T , respectively. The meaning of each column in these tables is as follows: “n” denotes the dimension of the problem; “|J0 |” stands for the number of elements in J0 ; “InIt” stands for the total number of inner iterations, namely, the total number of CG-iterations to
926
D.-H. Li et al. / Mathematical and Computer Modelling 48 (2008) 918–928
Table 1 Test results for Algorithm 1 with x0 = (1, . . . , 1)T n
InIt
OutIt
|J0 |
CPU
100 400 900 1600 2500 3600 4900 6400 8100 10 000 12 100 14 400 16 900 19 600 22 500
26 53 83 118 156 193 230 264 292 314 340 371 395 420 445
4 6 8 9 11 14 15 18 19 20 22 23 25 26 28
94 385 872 1551 2426 3497 4764 6222 7877 9730 11 775 14 014 16 451 19 084 21 906
0.063 0.172 0.343 0.703 1.313 2.234 3.688 5.515 7.938 11.187 14.594 19.766 24.938 30.672 37.375
Table 2 Test results for Algorithm 1 with x0 = (0, . . . , 0)T n
InIt
OutIt
|J0 |
CPU
100 400 900 1600 2500 3600 4900 6400 8100 10 000 12 100 14 400 16 900 19 600 22 500
20 40 63 93 125 157 187 213 237 256 285 305 323 338 364
4 6 8 9 11 14 15 17 18 20 22 23 25 26 28
94 385 872 1551 2426 3497 4764 6222 7877 9730 11 775 14 014 16 451 19 084 21 906
0.046 0.14 0.281 0.578 1.11 1.937 3.25 4.969 6.75 9.172 12.734 17 20.875 25.953 32.063
solve linear equations (3.6); “OutIt” means the total number of outer iterations and “CPU” denotes the computer time used (in second); We see from Tables 1 and 2 that starting from each initial point, the method terminates at the solution of the problem successfully. The proposed method can solve large-scale problems efficiently. 5. Discussion We have developed a conjugate gradient method for solving the linear complementarity problem with S-matrix. The method is an extension of the conjugate gradient method for solving systems of linear equations and minimization problems. The method reserves good properties of the conjugate gradient method for solving system of linear equations and minimization problems. An attractive property of the proposed method is that the computation cost of the method is no more than the cost of a conjugate gradient method for solving an n-dimensional system of linear equations. The proposed method depends on useful properties of S-matrix. These properties were used to identify the positive variables. It is interesting to develop a conjugate gradient method for solving the linear complementarity problem with general symmetrical and positive definite matrix. The difficulty to do this lies in the identification of the positive
D.-H. Li et al. / Mathematical and Computer Modelling 48 (2008) 918–928
927
variables. It is also interesting to develop preconditioned conjugate gradient method. This seems not trivial because Algorithm 1 implicitly implies the solution of a sequence of systems of linear equations. We leave it as a topic for further research. Acknowledgements The authors would like to thank the two anonymous referees for their very helpful comments. Appendix The Proof of Lemma 2.2. Since M is an M-matrix, its principal submatrices Mαα and Mββ are also M-matrices (see, e.g., in [13]), and then the Schur complement Mββ − Mβα (Mαα )−1 Mαβ is a Z -matrix. Noting that Iαα 0 Mαα Mαβ Iαα −(Mαα )−1 Mαβ Mβα Mββ −Mβα (Mαα )−1 Iββ 0 Iββ Mαα 0 = 0 Mββ − Mβα (Mαα )−1 Mαβ and Mββ − Mβα (Mαα )−1 Mαβ = Bββ , where
Bαα Bβα
Bαβ Bββ
=
Mαα Mβα
Mαβ Mββ
−1 ≥ 0,
we have that Mββ − Mβα (Mαα )−1 Mαβ is nonsingular and [Mββ − Mβα (Mαα )−1 Mαβ ]−1 ≥ 0. Hence, Mββ − Mβα (Mαα )−1 Mαβ is an M-matrix and then M¯ is an M-matrix by (2.5). Therefore, LCP(M, q) and ¯ q) LCP( M, ¯ exist unique solution respectively. So, it is sufficient to prove that the solution of LCP(M, q) is also the ¯ q). solution of LCP( M, ¯ ∗ Let x = x ≥ 0 be the solution of LCP(M, q). Since α ⊂ J (x ∗ ), by Lemma 2.1, we have m iT x − qi = 0,
i ∈ α,
i.e., Mαα xα + Mαβ xβ − qα = 0. Therefore, xα + (Mαα )−1 Mαβ xβ − (Mαα )−1 qα = 0.
(A.8)
And then (M x − q)β = Mβα xα + Mββ xβ − qβ = Mβα (Mαα )−1 (qα − Mαβ xβ ) + Mββ xβ − qβ = [Mββ − Mβα (Mαα )−1 Mαβ ]xβ − [qβ − Mβα (Mαα )−1 qα ]. Since x =
x∗
≥ 0 is the solution of LCP(M, q), by (A.8) and (A.9), we get
xα ≥ 0,
( M¯ x − q) ¯ α = xα + (Mαα )−1 Mαβ xβ − (Mαα )−1 qα = 0,
xβ ≥ 0,
( M¯ x − q) ¯ β = (M x − q)β ≥ 0,
and xβT ( M¯ x − q) ¯ β = 0,
¯ q). which implies x is the solution of LCP( M, ¯ The proof is then completed.
(A.9)
928
D.-H. Li et al. / Mathematical and Computer Modelling 48 (2008) 918–928
References [1] J.R. Birge, A. Gana, Computational complexity of Van der Heyden’s variable dimension algorithm and Dantzig–Cottle’s principal pivoting method for solving LCP’s, Mathematical Programming 26 (1983) 316–325. [2] R. Chandrasekaran, A special case of the complementarity pivot problem, Opsearch 7 (1970) 263–268. [3] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland Publishing Company, Amsterdam, 1978. [4] R.W. Cottle, J-S. Pang, R.E. Stone, The Linear Complementarity Problem, Academic Press, Inc., 1992. [5] C.W. Cryer, The eficient solution of linear complementarity problems for tradiagonal Minkowski matrices, ACM Transactions on Mathematical Software 9 (1983) 199–214. [6] F. Facchinei, A. Fischer, C. Kanzow, On the accurate identification of active constraints, SIAM Journal on Optimization 9 (1999) 14–32. [7] W.W. Hager, H. Zhang, A new active set algorithm for box constrained optimization, SIAM Journal on Optimization 17 (2007) 526–557. [8] H. Kremers, D. Talman, A new pivoting algorithm for the complementarity problem allowing for an arbitrary starting point, Mathematical Programming 63 (1994) 235–252. [9] D.H. Li, J.P. Zeng, Z.Z. Zhang, Gaussian pivoting method for solving linear complementarity problem, Applied Mathematics - A Journal of Chinese Universities 12 (1997) 419–426. [10] O.L. Mangasarian, Linear complementarity problems solvable by A single linear program, Mathematical Programming 10 (1976) 263–270. [11] K.G. Murty, Linear Complementarity, Linear and Nonlinear Programming, Heldermann, Berlin, 1988. [12] J. Nocedal, S. Wright, Numerical Optimization, Springer-Verlag, 2000. [13] J.M. Ortega, W.C. Rheinboldt, Iterative Solution of Nonlinear Equations in Several Variables, Academic Press, New York, 1970. [14] J. Rohn, A short proof of finiteness of Murty’s principal pivoting algorithm, Mathematical Programming 46 (1990) 255–256. [15] R. Saigal, Lemke’s algorithm and a special class of linear complementarity problems, Opsearch 8 (1971) 201–208. [16] M.J. Todd, Polynomial expected behavior of a pivoting algorithm for linear complementarity and linear programming problems, Mathematical Programming 35 (1986) 173–192. [17] L. Zhang, X.Y. Hu, On the direct method for linear complementarity problem, Mathematik Numerica Sinica 1 (1994) 59–64. [18] S.Z. Zhou, A direct method for the linear complementarity problem, Journal of Computational Mathematics 2 (1990) 178–182.