A generalization of the inexact parameterized Uzawa methods for saddle point problems

Applied Mathematics and Computation 206 (2008) 765–771

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

A generalization of the inexact parameterized Uzawa methods for saddle point problems Fang Chen a,b,*, Yao-Lin Jiang a a b

Department of Mathematical Sciences, Xi’an Jiaotong University, Xi’an 710049, PR China Department of Mathematics and Physics, Xi’an University of Post and Telecommunication, Xi’an 710121, PR China

a r t i c l e

i n f o

a b s t r a c t For large sparse saddle point problems, Bai and Wang recently studied a class of parameterized inexact Uzawa methods (see Z.-Z. Bai, Z.-Q. Wang, On paramaterized inexact Uzawa methods for generalized saddle point problems, Linear Algebra Appl. 428 (2008) 2900–2932). In this paper, we generalize these methods and propose a class of generalized inexact parameterized iterative schemes for solving the saddle point problems. We derive conditions for guaranteeing the convergence of these iterative methods. With different choices of the parameter matrices, the generalized iterative methods lead to a series of existing and new iterative methods including the classical Uzawa method, the inexact Uzawa method, the GSOR method and the GIAOR method. Ó 2008 Elsevier Inc. All rights reserved.

Keywords: Saddle point problems Iterative method Convergence

1. Introduction We consider solutions of the block 2  2 linear systems of the form

" AX ¼

A

BT

#  x

B

0

y

¼

  f ; g

ð1Þ

where A 2 Rnn is symmetric positive definite, B 2 Rmn ðm 6 nÞ is of full row rank, namely, rankðBÞ ¼ m x; f 2 Rn and y; g 2 Rm . We use BT to denote the transpose of the matrix B. This linear system is called saddle point problem. We know that the solution of (1) exists and is unique under the above mentioned assumptions; see [2–5]. Saddle point problems arise in many scientific computing and engineering applications such as constrained optimization, computational fluid dynamics, mixed finite element methods for solving elliptic partial differential equations and Stokes problems, constrained least-squares problems, structure analysis and so forth; see [15,13,12,3]. For this kind of problems, the best known and the oldest methods are the Uzawa and the preconditioned Uzawa methods. In recent years, a large variety of methods for solving this class of linear systems can be found in the literature. For example, Uzawa-type scheme [5], GSOR [2], HSS [3,6–8] and the Krylov subspace methods [9,10], etc. When the matrix blocks A and B are large and sparse, iterative methods become more attractive than direct methods for solving the saddle point problem (1), but direct methods play an important role in the form of preconditioners embedded in an iterative framework; see [2]. In this paper, we consider more generalized iterative methods by choosing different parameter matrices, and also give the corresponding convergence analysis. The remainder of the paper is organized as follows. In Section 2, we propose a generalization of the parameterized inexact Uzawa method for saddle point problems, and derive the condition for guaranteeing its convergence in Section 3. In Section 4, * Corresponding author. Address: Department of Mathematical Sciences, Xi’an Jiaotong University, Xian Ning West Road 28, Xi’an 710049, PR China. E-mail address: [email protected] (F. Chen). 0096-3003/$ - see front matter Ó 2008 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2008.09.041

766

F. Chen, Y.-L. Jiang / Applied Mathematics and Computation 206 (2008) 765–771

we use numerical examples to show the fast convergence of the generalized iterative method, which also show that the new method is an efficient and powerful tool for solving the large sparse saddle point problems. 2. Generalized iterative method The properties of the nonsymmetric matrix A are summarized in the following lemmas. Lemma 2.1. Let A be the coefficient matrix given in (1), where A is symmetric positive definite and B has full row rank. Let rðAÞ denote the spectrum of A. Then r A is nonsingular and detðAÞ > 0; s A is positive real, i.e., nT An > 0 holds for all n 2 Rnþm ðn – 0Þ; t A is positive stable, i.e., ReðkÞ > 0 holds for all k 2 rðAÞ. We consider the following matrix splitting:

"

A

BT

B

0

#

 ¼

"



A þ Q1

0

B þ Q 3

Q2



Q1

BT

Q3

Q2

# ;

where A þ Q 1 and Q 2 are nonsingular, and Q 3 2 Rmn is arbitrary. Then the generalized iterative method can be defined as



A þ Q1

0

B þ Q 3

Q2

"

xnþ1

#

ynþ1

"

¼

Q1

BT

Q3

Q2

#

xn



yn

þ

  f ; g

ð2Þ

or equivalently, it can be written as

(

xnþ1 ¼ xn þ ðA þ Q 1 Þ1 ðf  Axn  BT yn Þ; nþ1 ynþ1 ¼ yn þ Q 1 þ Q 3 xn þ gÞ: 2 ððB  Q 3 Þx

ð3Þ

Note that the iterative matrix of this iterative scheme is

W¼



A þ Q1

0

B þ Q 3

Q2

1 "

Q1 Q3

# BT : Q2

ð4Þ

Let qðWÞ denote the spectral radius of the iterative matrix W. Then the iterative scheme (3) converges when qðWÞ < 1; see [2,5–8]. Let k be an eigenvalue of W and ðuT ; vT ÞT be a corresponding eigenvector, where u 2 C n and v 2 C m . Then we have

W

    u u ¼k ; v v

or equivalently,

(

Q 1 u  BT v ¼ kðA þ Q 1 Þu; Q 3 u þ Q 2 v ¼ k½ðB þ Q 3 Þu þ Q 2 v:

ð5Þ

To prove the convergence of the iterative scheme (3), we need a few useful lemmas. Lemma 2.2. Assume that A is symmetric positive definite, B has full row rank, A þ Q 1 and Q 2 are nonsingular, and Q 3 2 Rmn is arbitrary. If k is an eigenvalue of the matrix W defined by (4), then k – 1. Proof. If k ¼ 1, then we have

(

BT v ¼ Au; Bu ¼ 0:

Note that these equations can be rewritten as

"

A

BT

#  u

B

0

v

¼ 0:

Clearly, the coefficient matrix is nonsingular. Hence, ðuT ; vT ÞT ¼ 0. This is a contradiction. So k – 1. h Lemma 2.3. Assume that A is symmetric positive definite, B has full row rank, A þ Q 1 and Q 2 are nonsingular, and Q 3 2 Rmn is arbitrary. If ðuT ; vT ÞT is an eigenvector of W with u 2 C n and v 2 C m , then u – 0.

F. Chen, Y.-L. Jiang / Applied Mathematics and Computation 206 (2008) 765–771

767

Proof. If u ¼ 0, then from (5) we have

(

BT v ¼ 0;

ð6Þ

ð1  kÞQ 2 v ¼ 0: Since k – 1 and Q 2 is nonsingular, by (6) we have v ¼ 0, which is a contradiction.

h

Lemma 2.4. Assume that A is symmetric positive definite, B has full row rank, Q 2 is nonsingular, Q 1 is symmetric positive semidefinite, and Q 3 2 Rmn is arbitrary. Let k be an eigenvalue of W and ðuT ; vT ÞT be a corresponding eigenvector, with u 2 C n and v 2 C m . If v ¼ 0, then 0 6 k < 1. Proof. If v ¼ 0, then from (5), we have



Q 1 u ¼ kðA þ Q 1 Þu;

ð7Þ

ð1  kÞQ 3 u ¼ kBu: Because u – 0, we can let

a¼

uH Q 1 u P 0; uH u

b¼

uH Au > 0: uH u

From the first equation in (7), we obtain

06k¼

a < 1:  aþb

Theorem 2.1. Assume that A is symmetric positive definite and B has full row rank. If Q 1 is symmetric positive semidefinite, Q 2 is symmetric positive definite, and Q 3 2 Rmn is such that BT Q 1 2 Q 3 is symmetric, Then the generalized iterative method (3) is convergent if and only if

c  4a  2b < 2s < 2b;

ð8Þ

where

a¼

uH Q 1 u P 0; uH u

b¼

uH Au > 0; uH u

c¼

uH BT Q 1 2 Bu P 0; uH u

s¼

uH BT Q 1 2 Q 3u uH u

and ðuT ; vT ÞT is an eigenvector of W, with u 2 C n and v 2 C m . Proof. Let k be an eigenvalue of W and ðuT ; vT ÞT be a corresponding eigenvector. From Lemmas 2.2 and 2.3, we have k – 1 and u – 0. So (5) can be rewritten as

(

Q 1 u  kðA þ Q 1 Þu ¼ BT v;

ð9Þ

1 v ¼ ðk1Þ Q 1 2 ½Q 3 u þ kðB  Q 3 Þu:

Eliminating t in (9), we get

½Q 1  kðA þ Q 1 Þu ¼

1 BT Q 1 2 ½Q 3 þ kðB  Q 3 Þu: ðk  1Þ

Multiplying both sides of this equality from left with

uH , uH u

we obtain

ðk  1Þ½a  kða þ bÞ ¼ s þ kðc  sÞ;

ð10Þ

where

a¼

uH Q 1 u P 0; uH u

b¼

uH Au > 0; uH u

c¼

uH BT Q 1 2 Bu P 0; uH u

s¼

uH BT Q 1 2 Q 3u : uH u

If c ¼ 0, then we only get Bu ¼ 0. When Bu ¼ 0, multiplying both sides of the first equation in (9) from left with uH , we have

06k¼

a < 1: aþb

When Bu – 0, we obtain c > 0. In this situation, (10) can be rewritten as

k2 

2a þ b þ s  c aþs kþ ¼ 0: aþb aþb

By [14,2,5], we know that a sufficient and necessary condition for the root of the real quadratic equation

x2 þ bx þ c ¼ 0

ð11Þ

768

F. Chen, Y.-L. Jiang / Applied Mathematics and Computation 206 (2008) 765–771

to satisfy jxj < 1 is that

jcj < 1 and jbj < 1 þ c:

ð12Þ

Hence, from (11) and (12) we know that jkj < 1 if and only if



ja þ sj < a þ b; j2a þ b þ s  cj < 2a þ b þ s:

Solving the above equations we obtain

c  4a  2b < 2s < 2b; where we have noticed that a P 0, b > 0 and c > 0.

h

Denote the smallest and the largest eigenvalues of A by kmin ðAÞ and kmax ðAÞ. Then we can get the following corollary. Corollary 2.1. Under the assumptions in Theorem 2.1, the generalized T 1 T 1 2BT Q 1 2 Q 3 þ 2A þ 4Q 1  B Q 2 B and A  B Q 2 Q 3 are positive definite.

iterative method (3)

is

convergent if

Proof. Based upon the proof of Theorem 2.1, we can get the following stronger convergence condition:

(

T 1 uH ð2BT Q 1 2 Q 3 þ 2A þ 4Q 1  B Q 2 BÞu > 0;

uH ðA  BT Q 1 2 Q 3 Þu > 0; T 1 T 1 which are implied in the conditions that both 2BT Q 1 h 2 Q 3 þ 2A þ 4Q 1  B Q 2 B and A  B Q 2 Q 3 are positive definite.

Corollary 2.2. Under the assumptions in Theorem 2.1, the generalized iterative method (3) is convergent if

(

kmax ðBT Q 1 2 Q 3 Þ < kmin ðAÞ; 1 k ðBT Q 1 2 BÞ 2 max

 2kmin ðQ 1 Þ  kmin ðAÞ < kmin ðBT Q 1 2 Q 3 Þ:

ð13Þ

Proof. From (8), we know

8 T 1 T 1 > < kmin ðB Q 2 Q 3 Þ 6 s 6 kmax ðB Q 2 Q 3 Þ; kmin ðAÞ 6 b; > :1 c  2a  b 6 12 kmax ðBT Q 1 2 BÞ  2kmin ðQ 1 Þ  kmin ðAÞ: 2 Hence the generalized iterative method is convergent if (13) holds. h If we specially choose Q 3 ¼ tB, then the generalized iterative method (2) becomes



A þ Q1 0 ðt  1ÞB Q 2

"

xnþ1

#

ynþ1

"

¼

Q1

BT

tB

Q2

#

   xn f þ : yn g

ð14Þ

Hence, the iterative method (14) is convergent when

c  4a  2b < 2tc < 2b: So a sufficient condition for the iterative method (14) to be convergent is T 1 T 1 kmax ðBT Q 1 2 BÞ  4kmin ðQ 1 Þ  2kmin ðAÞ < 2tkmin ðB Q 2 BÞ < 2tkmax ðB Q 2 BÞ < 2kmin ðAÞ:

3. Several algorithms Under the assumption of Theorem 2.1, the iterative method (3) leads to different iterative methods with different choices of the matrices Q 1 ; Q 2 and Q 3 . Case I (see [1]): Q 1 ¼ 0; Q 2 ¼ 1d I; Q 3 ¼ 0ðd > 0Þ. This gives the following Uzawa method. Algorithm 3.1. see [1]

(

xnþ1 ¼ xn þ A1 ðf  Axn  BT yn Þ; ynþ1 ¼ yn þ dðBxnþ1 þ gÞ:

Case II (see [11]): Q 1 is symmetric positive semidefinite, Q 2 is symmetric positive definite and Q 3 ¼ 0. This gives the following inexact Uzawa method.

769

F. Chen, Y.-L. Jiang / Applied Mathematics and Computation 206 (2008) 765–771

Algorithm 3.2. see [11]

(

xnþ1 ¼ xn þ ðA þ Q 1 Þ1 ðf  Axn  BT yn Þ; nþ1 þ gÞ: ynþ1 ¼ yn þ Q 1 2 ðBx

Case III (see [2]): Q 1 ¼ 1d A  A; Q 2 ¼ 1n Q ; Q 3 ¼ 0ðd > 0Þ and Q is symmetric positive definite. This gives the following generalized successive overrelaxation (GSOR) method. Algorithm 3.3. see [2]

(

xnþ1 ¼ ð1  dÞxn þ dA1 ðf  BT yn Þ; ynþ1 ¼ yn þ nQ 1 ðBxnþ1 þ gÞ:

In particular, when n ¼ d, the GSOR method becomes the SOR-like method. Case IV (see [2]): Q 1 ¼ 1d P  A; Q 2 ¼ 1n Q ; Q 3 ¼ B  nf Bðd > 0Þ, P  A; Q  BP1 BT , and P is symmetric positive definite. This gives the following special case of the GIAOR method. Algorithm 3.4. see [2]

(

xnþ1 ¼ xn þ dP1 ðf  Axn  BT yn Þ; ynþ1 ¼ yn þ nQ 1 ðBxn þ gÞ þ fQ 1 Bðxnþ1  xn Þ:

T Case V: Q 1 ¼ dA; Q 2 ¼ BQ 1 1 B and Q 3 ¼ tBðt > 0Þ. This gives the following new iterative method.

Algorithm 3.5

(

1 xnþ1 ¼ xn þ ð1þdÞ A1 ðf  Axn  BT yn Þ; nþ1 þ tBxn þ gÞ: ynþ1 ¼ yn þ Q 1 2 ðð1  tÞBx

4. Numerical examples In this section, we use an example from [7] to examine the numerical feasibility and effectiveness of Algorithm 3.5. All our tests are performed in MATLAB with machine precision 1016 , and terminated when the current residual satisfies kr k k2 =kr0 k2 < 105 , where rk is the residual at the kth iteration. We use the zero vector as the initial guess. In actual computations, we choose the right-hand-side vector ðf T ; g T ÞT such that the exact solution of the saddle point problem (1) is ðxT ; yT ÞT ¼ ð1; 1; . . . ; 1ÞT . We use IT and CPU to represent the number of iteration steps and the elapsed CPU time in seconds, respectively. In the numerical tables we denote by ERR ¼ krk k2 =kr0 k2 . Example 4.1 ([7,2]). Consider the saddle point problem (1), in which

 A¼

IT þT I

0

0

IT þT I

 ;

BT ¼



IF



FI

and

T¼

1 h

2

 tridiagð1; 2; 1Þ 2 Rpp ;

F¼

1  tridiagð1; 1; 0Þ 2 Rpp h

1 with  being the Kronecker product symbol and h ¼ pþ1 the discretization meshsize. Note that n ¼ 2p2 and m ¼ p2 .

Algorithm 3.5 is employed to solve Example 4.1. We list the computed results in Tables 1–5. The results for Algorithms 3.1 and 3.3 are listed in Tables 6–7 where the optimal parameter are employed. We observe from these tables that with the optimal choices of d and t, the convergence rate of Algorithm 3.5 is faster than the other algorithms. Evidently, we see that Algorithm 3.5 is much effective. The numerical results are in agreement with the theoretical one in Theorem 2.1. Table 1 d ¼ 0:4; t ¼ 0:01. n

m

IT

CPU

ERR

128 512 1152

64 256 576

18 18 18

0.063 1.703 11.344

8:9796  106 8:4478  106 8:2386  106

770

F. Chen, Y.-L. Jiang / Applied Mathematics and Computation 206 (2008) 765–771

Table 2 d ¼ 0:6; t ¼ 0:01. n

m

IT

CPU

ERR

128 512 1152

64 256 576

23 23 23

0.063 2.203 15.64

2:1113  106 2:1376  106 2:1679  106

Table 3 d ¼ 0:8; t ¼ 0:01. n

m

IT

CPU

ERR

128 512 1152

64 256 576

28 25 25

0.094 2.344 16.391

8:8231  106 9:9684  106 9:3251  106

n

m

IT

CPU

ERR

128 512 1152

64 256 576

29 28 28

0.094 2.703 17.922

9:9594  106 9:5681  106 9:4156  106

n

m

IT

CPU

ERR

128 512 1152

64 256 576

29 29 29

0.125 2.735 19.593

8:7966  107 4:7027  107 3:1251  107

n

m

IT

CPU

ERR

128 512 1152

64 256 576

24 34 42

0.063 2.203 20.407

8:6286  106 9:3683  106 9:5738  106

n

m

IT

CPU

ERR

128 512 1152

64 256 576

22 31 38

0.047 2.422 19.359

4:7596  107 6:1963  107 6:8247  107

Table 4 d ¼ 0:8; t ¼ 0:02.

Table 5 d ¼ 0:8; t ¼ 0:03.

Table 6 Algorithm 3.1

Table 7 Algorithm 3.3

References [1] K. Arrow, L. Hurwicz, H. Uzawa, Studies in Nonlinear Programming, Stanford University Press, Stanford, 1958. [2] Z.-Z. Bai, B.N. Parlett, Z.-Q. Wang, On generalized successive overrelaxation methods for augmented linear systems, Numer. Math. 102 (2005) 1–38. [3] Z.-Z. Bai, G.H. Golub, C.-K. Li, Optimal parameter in Hermitian and skew-Hermitian splitting method for certain two-by-two block matrices, SIAM J. Sci. Comput. 28 (2006) 583–603. [4] Z.-Z. Bai, J.-Y. Pan, M.K. Ng, New preconditioners for saddle point problems, Appl. Math. Comput. 172 (2006) 762–771. [5] Z.-Z. Bai, Z.-Q. Wang, On parameterized inexact Uzawa methods for generalized saddle point problems, Linear Algebra Appl. 428 (2008) 2900–2932. [6] Z.-Z. Bai, G.H. Golub, M.K. Ng, Hermitian and skew-Hermitian splitting methods for non-Hermitian positive definite linear systems, SIAM J. Matrix Anal. Appl. 24 (2003) 603–626. [7] Z.-Z. Bai, G.H. Golub, J.-Y. Pan, Preconditioned Hermitian and skew-Hermitian splitting methods for non-Hermitian positive semidefinite linear systems, Numer. Math. 98 (2004) 1–32.

F. Chen, Y.-L. Jiang / Applied Mathematics and Computation 206 (2008) 765–771

771

[8] Z.-Z. Bai, G.H. Golub, Accelerated Hermitian and skew-Hermitian splitting iteration methods for saddle point problems, IMA J. Numer. Anal. 27 (2007) 1–23. [9] Z.-Z. Bai, G.-Q. Li, Restrictively preconditioned conjugate gradient methods for systems of linear equations, IMA J. Numer. Anal. 23 (2003) 561–580. [10] Z.-Z. Bai, Structured preconditioners for nonsingular matrices of block two-by-two structures, Math. Comput. 75 (2006) 791–815. [11] J.H. Bramble, J.E. Pasciak, A.T. Vassilev, Analysis of the inexact Uzawa algorithm for saddle point problems, SIAM J. Numer. Anal. 34 (1997) 1072–1092. [12] H.C. Elman, D.J. Silvester, A.J. Wathen, Iterative methods for problems in computational fluid dynamics, in: Iterative Methods in Scientific Computing, Springer-Verlag, Singapore, 1997, pp. 271–327. [13] C.-J. Li, B.-J. Li, D.J. Evans, A generalized successive overrelaxation method for least squares problems, BIT 38 (1998) 347–356. [14] J.J.H. Miller, On the location of zeros of certain classes of polynomials with applications to numerical analysis, J. Inst. Math. Appl. 8 (1971) 397–406. [15] S. Wright, Stability of augmented system factorizations in interior-point methods, SIAM J. Matrix Anal. Appl. 18 (1997) 191–222.