- Email: [email protected]
On unsymmetric block overrelaxation-type methods for saddle point problems
Applied Mathematics and Computation 203 (2008) 660–671
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
On unsymmetric block overrelaxation-type methods for saddle point problems q Xiao-Fei Peng, Wen Li * School of Mathematical Sciences, South China Normal University, Guangzhou 510631, PR China
a r t i c l e
i n f o
a b s t r a c t The unsymmetric block overrelaxation-type (UBOR-type) method is proposed to attack saddle point problems in this paper. The convergence and the optimal parameters for the method are studied when the iteration parameters satisfy some relationship. Theoretical analyses show that the UBOR-type method has faster asymptotic convergence rate than the SOR-like method and its convergence rate can reach the same as that of the GSOR method at least. Numerical experiments support our theoretical results. Moveover, the numerical results further reveal that the new method can be much more effective than the GSOR method in terms of iteration steps. Ó 2008 Elsevier Inc. All rights reserved.
Keywords: The saddle point problem The UBOR-type method Convergence Optimal parameters
1. Introduction Consider the saddle point problem
A
B T
B
0
x y
¼
b q
;
or
eu ~ ¼ ~f ; A
ð1:1Þ
where A 2 Rmm is a symmetric positive definite matrix, B 2 Rmn is a matrix of full column rank with m > n and b 2 Rm , q 2 Rn are two given vectors. The saddle point problem appears frequently in many scientific and engineering computing [1–3]. A number of solution methods for solving the linear system (1.1) can be found in the literature, such as Uzawa-type schemes [4–7], splitting methods [8–14], preconditioned Krylov subspace methods [12,15–19]. Recently, Golub et al. proposed the SOR-like method and considered the optimum choice for the iterative parameter by using some nonsingular preconditioning e [9]. Li and Shao et al. generalized the relaxed splitting method to the matrix Q instead of the null block in coefficient matrix A GAOR method to improve the convergence rate [20,21]. However, they do not achieve the optimal parameters, which affects the efficiency of the method. By using a parameter matrix instead of a single parameter, Bai et al. presented the GSOR method and accelerated the SOR-like method without extra cost per iteration step [22]. Additionally, by using the techniques of vector extrapolation, matrix relaxation and inexact iteration, the GSOR method was extended to the generalized inexact accelerated overrelaxation (GIAOR) method in [22]. To further improve the convergence rate, we present an unsymmetric block overrelaxation-type (UBOR-type) method with three parameters s; x1 ; x2 in this paper. In Section 2, the new splitting method is proposed and its convergence properties are studied. In Section 3, the optimal parameters and corresponding spectral radius are determined when the param-
q
The work was supported by National Natural Science Foundation (Grant No. 10671077) and Guangdong provincial Natural Science Foundation (Grant No. 06025061). * Corresponding author. E-mail addresses: [email protected] (X.-F. Peng), [email protected] (W. Li).
0096-3003/$ - see front matter Ó 2008 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2008.05.014
661
X.-F. Peng, W. Li / Applied Mathematics and Computation 203 (2008) 660–671
eters satisfy the relationship s ¼ x1 þ x2 x1 x2 . Numerical experiments are presented to support our theoretical results in Section 4. Finally, we draw our conclusions in Section 5. 2. The UBOR-type method Let Q 2 Rnn be nonsingular and symmetric, and
e¼ A
A
B
BT
0
¼ D L U;
where
D¼
A 0 ; 0 Q
L¼
0 ; 0
0 BT
U¼
0 B : 0 Q
Then the UBOR-type iteration scheme for solving the augmented linear system (1.1) is given by
xðkþ1Þ
! ¼ Gs;x1 ;x2
yðkþ1Þ
xðkÞ
! þ M 1 s;x1 ;x2
yðkÞ
b ; q
and the iteration matrix Gs;x1 ;x2 of the UBOR-type method can be written as
Gs;x1 ;x2 ¼ M1 s;x1 ;x2 N s;x1 ;x2 ; where
M s;x1 ;x2 ¼
1
s
ðD x1 LÞD1 ðD x2 UÞ ¼
1
A T
s x1 B
x2 B ð1 x2 ÞQ x1 x2 BT A1 B
and
e¼ Ns;x1 ;x2 ¼ M s;x1 ;x2 A ¼
1
s
ð1 sÞA s ðs x1 ÞBT
1
½ð1 sÞD þ ðs x1 ÞL þ ðs x2 ÞU þ x1 x2 LD1 U ! ðx2 sÞB : ð1 x2 ÞQ x1 x2 BT A1 B
Using the block-inversion formula
A11
A12
A21
A22
1 ¼
1 1 1 ½I þ A1 11 A12 ðA22 A21 A11 A12 Þ A21 A11
1 1 A1 11 A12 ðA22 A21 A11 A12 Þ
1 1 ðA22 A21 A1 11 A12 Þ A21 A11
1 ðA22 A21 A1 11 A12 Þ
! ;
we obtain
0 ðsMs;x1 ;x2 Þ1 ¼ @
1 x2 A1 BQ 1 BT A1 A1 x 1x2
x1
1x2
Q 1 BT A1
1xx2 2 A1 BQ 1 1 1x2
Q 1
1 A:
Hence 1 Gs;x1 ;x2 ¼ M1 s;x ;x N s;x1 ;x2 ¼ ðsM s;x1 ;x2 Þ ðsN s;x1 ;x2 Þ 0 1 2 1 1 1 T 1 T 1 x1 Þx2 s 1 1 x2 s sA1 B þ x1 ð1 sÞI ð11 x2 A BQ B x2 A BQ B A B A @ ; ¼ 1 T 1 ð1x1 Þs 1 T x1 s Q B I 1 1x2 x2 Q B A B
where I is the identity matrix, parameters s; x1 ; x2 are real numbers with
ð2:1Þ
s 6¼ 0; x2 6¼ 1.
Algorithm 2.1 (The UBOR-type method) 1. Given initial vectors xð0Þ ; yð0Þ , the relaxation parameters s; x1 ; x2 , and the preconditioning matrix Q, which is an approximate matrix of the Schur complement matrix BT A1 B. 2. For k ¼ 0; 1; . . . ; until convergence, do 1 T 1 1 ðkÞ T ðkÞ sx1 s yðkþ1Þ ¼ yðkÞ þ 1 x2 Q B A ðb By Þ þ 1x2 Q ½ð1 x1 ÞB x q;
xðkþ1Þ ¼ ð1 sÞxðkÞ þ sA1 ðb ByðkÞ Þ x2 A1 Bðyðkþ1Þ yðkÞ Þ: It is worth mentioning that we have computed yðkþ1Þ before computing xðkþ1Þ . Obviously, when x2 ¼ 0, the above algorithm reduces to the GAOR method [20]:
662
X.-F. Peng, W. Li / Applied Mathematics and Computation 203 (2008) 660–671
xðkþ1Þ ¼ ð1 sÞxðkÞ þ sA1 ðb ByðkÞ Þ; yðkþ1Þ ¼ yðkÞ Q 1 fsq þ BT ½ðx1 sÞxðkÞ x1 xðkþ1Þ g: 1s When three parameters satisfy the relationship s ¼ x1 þ x2 x1 x2 , i.e., x1 ¼ 1 1 x2 , the UBOR-type method reduces to the following two-parameter UBOR-type (TP-UBOR-type) method:
sðsx2 Þ 1 T 1 yðkþ1Þ ¼ yðkÞ þ ð1 Q B A ðb ByðkÞ Þ þ 1sx2 Q 1 x Þ2
2
1s 1x2
BT xðkÞ q ;
xðkþ1Þ ¼ ð1 sÞxðkÞ þ sA1 ðb ByðkÞ Þ x2 A1 Bðyðkþ1Þ yðkÞ Þ: e namely, Remark 2.1. The TP-UBOR-type method is a fixed point iteration method based on the splitting of A, e ¼ Ms;x N s;x , where A 2 2
M s;x2 ¼
!
x2 B
A
1
ðsx2 Þ T 1 sx2 T ð1 x2 ÞQ x21 s 1 x2 B x2 B A B
and
Ns;x2 ¼
ð1 sÞA
1
!
ðx2 sÞB
ðsx2 Þ T 1 2 ð1sÞ T s x1 ð1 x2 ÞQ x21 x2 B x2 B A B
:
Specially, if x2 ¼ 0, the TP-UBOR-type algorithm reduces to the SOR-like method [9]:
xðkþ1Þ ¼ ð1 sÞxðkÞ þ sA1 ðb ByðkÞ Þ; yðkþ1Þ ¼ yðkÞ þ sQ 1 ðBT xðkþ1Þ qÞ: Remark 2.2. It is noted that the TP-UBOR-type method and GSOR method in [22] have different algorithmic structures based on different splittings of the coefficient matrix, but they can both become the SOR-like method under the different restrictions on the iteration parameters. Generally speaking, the computational cost of the former can be more than that of the latter at each step. To improve the running speed, we consider the following practical Algorithm. Algorithm 2.2 (The practical UBOR-type method) 1. Given initial vectors xð0Þ ; yð0Þ , the relaxation parameters s; x1 ; x2 , and the preconditioning matrix Q, which is an approximate matrix of the Schur complement matrix BT A1 B. 2. Let t1 ¼ 1sx2 ; t 0 ¼ t1 x1 ; t 2 ¼ 1 s; t3 ¼ 1 x1 . 3. Let T 1 ¼ Q 1 BT ; T 2 ¼ A1 b; T 3 ¼ A1 B; T 4 ¼ Q 1 q. 4. For k ¼ 0; 1; . . . ; until convergence, do
T ¼ T 2 T 3 yðkÞ ; yðkþ1Þ ¼ yðkÞ þ t0 T 1 T þ t 1 ½t 3 T 1 xðkÞ T 4 ; xðkþ1Þ ¼ t2 xðkÞ þ ðs x2 ÞT þ x2 ðT 2 T 3 yðkþ1Þ Þ: Next, we shall study the convergence of the UBOR-type method for soving the linear system (1.1). By rðÞ and qðÞ we denote the spectrum and spectral radius of a matrix, respectively. Let k 2 rðGs;x1 ;x2 Þ and l 2 rðQ 1 BT A1 BÞ; lmin ¼ minflg and q ¼ qðQ 1 BT A1 BÞ. Suppose that ðxT ; yT ÞT is the eigenvector corresponding to k. Then
Gs;x1 ;x2
x y
x ¼k ; y
i.e.,
ð1 sÞA
ðs x2 ÞB
ðs x1 ÞBT
ð1 x2 ÞQ x1 x2 BT A1 B
! x y
¼k
A x1 BT
x2 B ð1 x2 ÞQ x1 x2 BT A1 B
x ; y
ð2:2Þ
from which one may deduce that
ð1 s kÞx ¼ ðkx2 þ s x2 ÞA1 By
ð2:3Þ
ð1 kÞ½x1 x2 BT A1 By ð1 x2 ÞQy ¼ ðkx1 þ s x1 ÞBT x:
ð2:4Þ
and
X.-F. Peng, W. Li / Applied Mathematics and Computation 203 (2008) 660–671
663
Lemma 2.1 [23]. Both roots of the real quadratic equation x2 bx þ c ¼ 0 are less than one in modulus if and only if jcj < 1 and jbj < 1 þ c: Lemma 2.2. Assume that 1 s 6¼ 0. Then (1) when x1 ¼ 1; k ¼ 1 s is an eigenvalue with m multiplicities at least; (2) when x1 6¼ 1 and m > n, k ¼ 1 s is an eigenvalue with m n multiplicities at least. Proof (1) When x1 ¼ 1, by (2.1) we obtain
0 Gs;x1 ;x2 ¼ @
ð1 sÞI
1 T 1 x2 s 1 sA1 B þ 1 x2 A BQ B A B
O
I 1sx2 Q 1 BT A1 B
1 A;
which implies the conclusion (1). (2) From (2.3) and (2.4), an explicit calculation reveals that k ¼ 1 s is an eigenvalue of Gs;x1 ;x2 with m n multiplicities at least and the corresponding eigenvector is ðxT ; 0ÞT , where ð1 x1 ÞBT x ¼ 0: In terms of the assumptions on x1 and B, it holds that x 2 NðBT Þ and x 6¼ 0. h Theorem 2.1. Let A 2 Rmm and Q 2 Rnn be symmetric positive definite, and B 2 Rmn be of full column rank. Then the UBOR-type method is convergent if and only if
0 < s < 2; x2 < 1
qs2
; 4 s x2 1 s 2x2 2 s < x1 < : þ sq 1 x2 q 2ð1 x2 Þ 1
1
Proof. Since Q 1 BT A1 B is similar to Q 2 BT A1 BQ 2 ; they have same eigenvalues. Moreover, the latter is congruent to the positive definite matrix BT A1 B: It follows from Sylvester’s law of inertia that l > 0 [24]. By Lemma 2.2, k ¼ 1 s 6¼ 0 is an eigenvalue of Gs;x1 ;x2 : If k 6¼ 1 s, it follows from (2.3) and (2.4) that
k2 ðsl slx1 þ 2 s
sl 1 x2
Þk þ sl slx1 þ 1 s
slð1 sÞ ¼ 0; 1 x2
ð2:5Þ
or equivalently,
k2 ðc þ 1
s2 l Þk þ c ¼ 0; 1 x2
where
c ¼ sl slx1 þ 1 s
slð1 sÞ : 1 x2
By Lemma 2.1, jkj < 1 if and only if
j1 sj < 1; jcj < 1;
2
s l jc þ 1 1 x2 j < 1 þ c:
It follows that
0 < s < 2;
x2 < 1;
and
s2 l 2ð1 x2 Þ
< c þ 1 < 2;
s2 l 2ð1 x2 Þ
< 2:
In terms of (2.2), we have
s x2 1 s x2 1 6 < x1 ; 1 x2 l 1 x2 q s 2x2 2 s s 2x2 2 s x1 < 6 : þ þ sq sl 2ð1 x2 Þ 2ð1 x2 Þ s2 l 6 s2 q < 4ð1 x2 Þ;
ð2:6Þ
664
X.-F. Peng, W. Li / Applied Mathematics and Computation 203 (2008) 660–671
Therefore, we obtain the following necessary and sufficient conditions such that the UBOR-type method is convergent
0 < s < 2; x2 < 1
qs2
; 4 s x2 1 s 2x2 2 s < x1 < : þ sq 1 x2 q 2ð1 x2 Þ
Let x2 ¼ 0; then the following result follows immediately from Theorem 2.1. Corollary 2.1 [20]. Under the same assumptions as in Theorem 2.1, the GAOR method is convergent if and only if
1 s 2s s < x1 < þ : sq q 2
2 0 < s < minf2; pffiffiffiffig;
q
Remark 2.3. Actually, the above GAOR method is a special GIAOR method in [22] since the latter may reduce to the former when we take P ¼ A, D ¼ cI and X ¼ xI. If
s ¼ x1 þ x2 x1 x2 , by Theorem 2.1 then we have:
Corollary 2.2. Under the same assumptions as in Theorem 2.1, the TP-UBOR-type method is convergent if and only if
0 < s < 2;
x2 < 1
q s2 2ð2 sÞ
:
Corollary 2.3 [9]. Under the same assumptions as in Theorem 2.1, the SOR-like method is convergent if and only if
4 0 < s < pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : 4q þ 1 þ 1 3. The optimal iteration parameters By the proof of Theorem 2.1, k ¼ 1 s or
k1;2
1 ¼ 2
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi)
(
sl l 2 4l ; sl slx1 þ 2 s s ðlx1 þ 1 l þ Þ 1 x2 1 x2 1 x2
ð3:1Þ
where k1 and k2 correspond to taking the ‘+’ and ‘’ signs in (3.1), respectively. Then, we have
qðGs;x1 ;x2 Þ ¼ max fj1 sj; jk1 j; jk2 jg: l According to Theorem 2.1, the parameters of the convergent UBOR-type method satisfy the following conditions:
0 < s < 2; x2 < 1
qs2
; 4 s x2 1 s 2x2 2 s < x1 < : þ sq 1 x2 q 2ð1 x2 Þ
To achieve the absolute values of the nonzero eigenvalues of the iteration matrix Gs;x1 ;x2 , we consider the following three cases: (1)
sl slx1 þ 2 s
(2)
sl slx1 þ 2 s
(3) ðlx1 þ 1 l þ
sl 1 x2
sl 1 x2
l 1 x2
Þ2
P 0;
ðlx1 þ 1 l þ
6 0;
ðlx1 þ 1 l þ
l 1 x2
l 1 x2
Þ2
4l P 0; 1 x2
Þ2
4l P 0; 1 x2
4l 6 0: 1 x2
Accordingly, we define the following functions, whose variables s; x1 ; x2 conform to the restrictions of the above three cases, respectively.
X.-F. Peng, W. Li / Applied Mathematics and Computation 203 (2008) 660–671
1 f1 ðs; x1 ; x2 ; lÞ ¼ 2 for
665
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi)
(
sl l 2 4l sl slx1 þ 2 s þ s ðlx1 þ 1 l þ Þ ; 1 x2 1 x2 1 x2
2 1 x2 1 2 x2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 6 x1 6 ð 1Þ ; l 1 x l s 1 x2 lð1 x2 Þ 2
s6
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 x2 pffiffiffiffi ;
l
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi) ( 1 sl l 2 4l f2 ðs; x1 ; x2 ; lÞ ¼ slx1 sl þ s 2 þ þ s ðlx1 þ 1 l þ Þ ; 2 1 x2 1 x2 1 x2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 x2 1 2 x2 1 ; s 6 pffiffiffiffi or for x1 P l s 1 x2 l pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 x2 2 1 x2 x1 P pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; s P pffiffiffiffi ; l lð1 x2 Þ l 1 x2 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi slð1 sÞ ; f3 ðs; x1 ; x2 ; lÞ ¼ sl slx1 þ 1 s 1 x2
2 1 x2 for x1 6 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : lð1 x2 Þ l 1 x2
Obviously,
qðGs;x1 ;x2 Þ ¼ j1 sj; or
8 2 ffi 1 x2 ; max f3 ðs; x1 ; x2 ; lÞ; for x1 6 pffiffiffiffiffiffiffiffiffiffiffiffiffi > > lð1x2 Þ l 1x2 l > > > > > 2 > ffi 1 x2 6 x1 6 l1 ð2s 1Þ 1xx2 ; max f1 ðs; x1 ; x2 ; lÞ; for pffiffiffiffiffiffiffiffiffiffiffiffiffi > > 2 lð1x2 Þ l 1x2 l > > > pffiffiffiffiffiffiffiffiffi > < 1x2 s 6 pffiffilffi ; qðGs;x1 ;x2 Þ ¼ > pffiffiffiffiffiffiffiffiffi > > 1x > > max f2 ðs; x1 ; x2 ; lÞ; for x1 P l1 ð2s 1Þ 1xx2 2 ; s 6 pffiffilffi 2 > > l > > > pffiffiffiffiffiffiffiffiffi > > > 1x2 2 > ffi l1 1xx2 ; s P p ffiffilffi : : or x1 P pffiffiffiffiffiffiffiffiffiffiffiffiffi 2
ð3:2Þ
lð1x2 Þ
Specially, we shall determine the optimal iteration parameters of the TP-UBOR-type method. Naturally, the notations
qðGs;x2 Þ and fj ðs; x2 ; lÞ in the TP-UBOR-type method can be used instead of qðGs;x1 ;x2 Þ and fj ðs; x1 ; x2 ; lÞ (j = 1, 2, 3) in the UBOR-type method, respectively. Theorem 3.1. Let A 2 Rmm and Q 2 Rnn be symmetric positive definite, B 2 Rmn be of full column rank. Then the spectral radius of the iteration matrix Gs;x2 of the TP-UBOR-type method is given as follows: (1) for 0 < s < 1;
8 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi s2 lmin slmin 2 4lmin 1 > > > 2 ½2 s 1x2 þ s ð1 þ 1x2 Þ 1x2 ; for > < pffiffiffiffiffiffiffiffiffiffiffiffi qðGs;x2 Þ ¼ for 1 s; > > qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > > : 1 ½s 2 þ s2 q þ s ð1 þ sq Þ2 4q ; for 1x2 1x2 1x2 2
x2 < x2ð1Þ ðsÞ; ð2Þ xð1Þ 2 ðsÞ 6 x2 6 x2 ðsÞ;
ð3:3Þ
2
s q xð2Þ 2 ðsÞ < x2 < 1 2ð2sÞ ;
where
xð1Þ 2 ðsÞ ¼ 1
s2 lmin
pffiffiffiffiffiffiffiffiffiffiffiffi ; ð1 1 sÞ2
x2ð2Þ ðsÞ ¼ 1
s2 q
pffiffiffiffiffiffiffiffiffiffiffiffi : ð1 þ 1 sÞ2
(2) for 1 6 s < 2; we have
8 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 slmin 2 4lmin > < 12 ½2 s s1lxmin þ s ð1 þ 1 x2 Þ 1x2 ; for 2 qðGs;x2 Þ ¼ q ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi > : 1 ½s 2 þ s2 q þ s ð1 þ sq Þ2 4q ; for 1x2 1x2 1x2 2 where
xð3Þ 2 ðsÞ ¼ 1
s2 ðq þ lmin Þ : 2ð2 sÞ
x2 < x2ð3Þ ðsÞ; 2
s q xð3Þ 2 ðsÞ 6 x2 < 1 2ð2sÞ ;
ð3:4Þ
666
X.-F. Peng, W. Li / Applied Mathematics and Computation 203 (2008) 660–671
Furthermore, the optimal parameters are given by
pffiffiffiffiffiffiffiffiffiffiffiffiffi
4 qlmin 4qlmin sb ¼ pffiffiffiffi pffiffiffiffiffiffiffiffiffi ffi 2 ; x2b ¼ 1 pffiffiffiffi pffiffiffiffiffiffiffiffiffiffi 2 ; ð q þ lmin Þ ð q þ lmin Þ and the corresponding optimal spectral radius of the iteration matrix Gs;x2 is
pffiffiffiffi pffiffiffiffiffiffiffiffiffiffi q lmin min qðGs;x2 Þ ¼ qðGsb ;x2b Þ ¼ pffiffiffiffi pffiffiffiffiffiffiffiffiffi ffi: s;x 2 q þ lmin
Proof. By Corollary 2.2, it need only to consider that parameters s and x2 satisfy
0 < s < 2;
x2 < 1
q s2 2ð2 sÞ
:
According to the expressions of the three functions fj ðs; x1 ; x2 ; lÞ ðj ¼ 1; 2; 3Þ we derive the following expressions of the functions fj ðs; x2 ; lÞ ðj ¼ 1; 2; 3Þ in the TP-UBOR-type method as follows:
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffio n s2 l sl 2 4l f1 ðs; x2 ; lÞ ¼ 12 2 s 1 x2 þ s ð1 þ 1x2 Þ 1x2 ; 2
l ; for 0 < s < 1; x2 < 1 ð1ps ffiffiffiffiffiffi 1sÞ2 2
s l or 1 6 s < 2; x2 < 1 2 s; qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffio n 2l s f2 ðs; x2 ; lÞ ¼ 12 s 2 þ 1x2 þ s ð1 þ 1slx2 Þ2 14lx2 ; 2
2
l s q < x2 < 1 2ð2 for 0 < s < 1; 1 ð1þps ffiffiffiffiffiffi sÞ ; 1sÞ2 2
2
s l s q or 1 6 s < 2; 1 2 s 6 x2 < 1 2ð2sÞ ; pffiffiffiffiffiffiffiffiffiffiffiffi f3 ðs; x2 ; lÞ ¼ 1 s; 2
2
l l 6 x2 6 1 ð1þps ffiffiffiffiffiffi : for 0 < s < 1; 1 ð1ps ffiffiffiffiffiffi 1sÞ2 1sÞ2
Under the above restrictions on s; x2 , it holds that f3 ðs; x2 ; lÞ ¼ fj ðs; x2 ; lÞ ðj ¼ 1; 2; 3Þ to achieve qðGs;x2 Þ. A simple computation reveals that
8 > > <
pffiffiffiffiffiffiffiffiffiffiffiffi 1 s P j1 sj: Hence, we just consider the functions
9 > > = of1 ðs; x2 ; lÞ s ¼ s þ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi < 0; > 2 ol 2ð1 x2 Þ > > sl 4l > : ; 1 þ 1x2 1x2 8 9 > > > > s2 l < = s 2 þ 1x2 of2 ðs; x2 ; lÞ s ¼ s þ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > 0; 2 > ol 2ð1 x2 Þ > > > : 1 þ sl 4l ; s2 l 2 s 1 x2
1x2
ð3:5Þ
1x2
of3 ðs; x2 ; lÞ ¼ 0; ol and
8 > > <
9 > > =
s2 l 2 s 1 x2
of1 ðs; x2 ; lÞ sl ffi < 0; ¼ s þ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 > o x2 2ð1 x2 Þ2 > > > : 1 þ sl 4l ; 8 > > <
1x2
1x2
9 > > =
2
s 2 þ 1s xl2
of2 ðs; x2 ; lÞ sl ffi > 0; ¼ s þ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 > o x2 2ð1 x2 Þ2 > > > : 1 þ sl 4l ; 1x2
1x2
of3 ðs; x2 ; lÞ ¼ 0: o x2 When 0 < s < 1; or 1 6 s < 2; then (3.3) or (3.4) follows immediately from (3.2) and (3.5). Now, we consider the optimal relaxation parameters in the following two cases. (1) For the case 0 < s < 1; it follows from (3.3) and (3.6) that
ð3:6Þ
667
X.-F. Peng, W. Li / Applied Mathematics and Computation 203 (2008) 660–671
8 ð1Þ ð1Þ > min f1 ðs; x2 ; lmin Þ ¼ f1 ðs; x2 ðsÞ; lmin Þ; x2 6 x2 ðsÞ; > > > x2 < pffiffiffiffiffiffiffiffiffiffiffiffi ð2Þ min qðGs;x2 Þ ¼ xð1Þ 1 s; 2 ðsÞ 6 x2 6 x2 ðsÞ; x2 > > > ð2Þ ð2Þ s2 q > : min f2 ðs; x2 ; qÞ ¼ f2 ðs; x2 ðsÞ; qÞ; x2 ðsÞ 6 x2 < 1 2ð2 sÞ :
ð3:7Þ
x2
For any fixed s, it is easy to verify that ð1Þ
ð2Þ
f1 ðs; x2 ðsÞ; lmin Þ ¼ f2 ðs; x2 ðsÞ; qÞ ¼
pffiffiffiffiffiffiffiffiffiffiffiffi 1 s: ð1Þ
Evidently, minx2 qðGs;x2 Þ attains the minimum when the parameter s attains the maximum. Here, it noted that x2 ðsÞ and ðiÞ xð2Þ 2 ðsÞ are the functions of s. Next, we investigate further the monotonicity of x2 ðsÞ; i ¼ 1; 2. Because
pffiffiffiffiffiffiffiffiffiffiffiffi ð1Þ ox2 ðsÞ lmin ð1 þ 1 sÞ pffiffiffiffiffiffiffiffiffiffiffiffi > 0; ¼ os 1s
pffiffiffiffiffiffiffiffiffiffiffiffi ð2Þ ox2 ðsÞ qð1 1 sÞ pffiffiffiffiffiffiffiffiffiffiffiffi < 0; ¼ os 1s
ð2Þ xð1Þ 2 ðsÞ is increasing and x2 ðsÞ is decreasing when s is increasing. To ensure that s both attains the maximum and satisfies ð1Þ
ð2Þ
ð1Þ
ð2Þ
the inequality x2 ðsÞ 6 x2 ðsÞ, the s such that qðGs;x2 Þ attains the minimum must satisfy the relation x2 ðsÞ ¼ x2 ðsÞ; namely,
1
s2 lmin
s2 q
pffiffiffiffiffiffiffiffiffiffiffiffi ¼ 1 pffiffiffiffiffiffiffiffiffiffiffiffi : ð1 1 sÞ2 ð1 þ 1 sÞ2
It then follows that ð2Þ xð1Þ qðGs;x2 Þ 2 ðsÞ ¼ x2 ðsÞ ¼ arg min x 2
for the fixed s. An simple calculation shows that the optimal relaxation parameters are
pffiffiffiffiffiffiffiffiffiffiffiffiffi
4 qlmin 4qlmin sbð1Þ ¼ pffiffiffiffi pffiffiffiffiffiffiffiffiffi ffi ; x2bð1Þ ¼ 1 pffiffiffiffi pffiffiffiffiffiffiffiffiffiffi 2 ; ð q þ lmin Þ ð q þ lmin Þ2 and the corresponding optimal spectral radius is
pffiffiffiffi
pffiffiffiffiffiffiffiffiffiffi
q lmin qðGðsbð1Þ ; x2bð1Þ ÞÞ ¼ pffiffiffiffi pffiffiffiffiffiffiffiffiffi ffi: q þ lmin (2) For the case 1 6 s < 2; it follows from (3.4) and (3.6) that ð3Þ
ð3Þ
min qðGs;x2 Þ ¼ f1 ðs; x2 ðsÞ; lmin Þ ¼ f2 ðs; x2 ðsÞ; qÞ x2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 ¼ ðq lmin Þð2 sÞ þ s2 ðq þ lmin Þ2 4qlmin ð2 sÞ2 : 2ðq þ lmin Þ Let yðsÞ ¼ 2ðq þ lmin Þ minx2 qðGðs; x2 ÞÞ: Because
oy ðq lmin Þ2 s þ 8qlmin P 0; ¼ ðq lmin Þ þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi os s2 ðq þ lmin Þ2 4qlmin ð2 sÞ2 the optimal relaxation parameters such that qðGðs; x2 ÞÞ attains the minimum are
sbð2Þ ¼ 1; x2bð2Þ ¼ xð3Þ 2 ¼ 1
q þ lmin 2
;
and the corresponding optimal spectral radius is
qðGsbð2Þ ;x2bð2Þ Þ ¼
q lmin : q þ lmin
Comparing the above two cases, the desired results follow. h Remark 3.1. It is noted that the sb of the optimal TP-UBOR-type method is just equal to the xb of the optimal GSOR method in [22] and their corresponding methods have the same asymptotic convergence factor. In other words, the UBOR-type method is equally effective to the GSOR method at least in theory since the TP-UBOR-type method is a special UBOR-type method.
668
X.-F. Peng, W. Li / Applied Mathematics and Computation 203 (2008) 660–671
Remark 3.2. Generally speaking, the optimal parameters sb ; x2b of the TP-UBOR-type method and the corresponding x1b are not optimal in the UBOR-type method, where the parameter x1b is determined by the relationship sb ¼ x1b þ x2b x1b x2b . By Theorem 3.1,
x1b ¼ 1
pffiffiffiffi
pffiffiffiffiffiffiffiffiffiffi
q lmin 2 : 4qlmin
It is very difficult to determine the optimal relaxation parameters of the UBOR-type method without the above given relationship on three parameters. However, as we shall show in next section on numerical experiments, the faster convergence rate of the UBOR-type method can be achieved for other cases different from those in Theorem 3.1.
4. Numerical experiments In this section, we apply the UBOR-type method to solve those problems taken from [7] and [22], and a set of Stokes saddle point problems generated using the IFISS software written by Howard Elman et al. Problem 1. Consider the linear system (1.1) with the coefficient sub-matrices A; B deriving from [7],
8 > < i þ 1; i ¼ j; A ¼ 1; ji jj ¼ 1; > : 0; other;
B¼
j;
i ¼ j þ m n;
0; other:
The right vectors are given by
b ¼ ð1; 1; . . . ; 1ÞT 2 Rm ;
q ¼ ð1; 1; . . . ; 1ÞT 2 Rn :
Problem 2. Consider the linear system with the coefficient sub-matrices A; B in [22],
A¼
IT þT I
0
0
IT þT I
;
B¼
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. The right vectors are defined as
b ¼ ð1; 1; . . . ; 1ÞT 2 Rm ;
q ¼ ð0; 0; . . . ; 0ÞT 2 Rn ; m ¼ 2p2 ; n ¼ p2 :
Problem 3. We generate the regularised-lid driven cavity problems with IFISS software using the 16 16; 32 32; 64 64 grid, Q 2 Q 1 discretisation. Since the matrix B generated by the package is actually singular, we drop the first two rows for the Stokes 1 problem and the first rows of B for the Stokes 2–3 problems to obtain the full rank matrices. The dimensions of the problems and number of nonzeros in the associated matrices A and B are given in Table 1. All experiments are performed in MATLAB with the machine precision 1016 and all numerical results are listed in Tables 2–7. To save the storage space, we just list the optimal parameters and spectral radius with six number of decimal place in 6 kr k k following Tables. We use a zero initial guess and the stopping criterion kr , where r k ; r 0 denote the kth residual vector 0 k 6 10 and the initial one, respectively. IT and CPU represent the number of iteration steps and the elasped CPU time in seconds, respectively. The preconditioning matrix Q is chosen as BT B for Problem 1 and is given as follows for Problems 2 and 3,
Case 1 : Q ¼ diagðBT A1 BÞ;
Case 2 : Q ¼ tridiagðBT A1 BÞ:
xb ; sb and x2b indicate the corresponding optimal parameters in the SOR-like, GSOR and UBOR-type methods. x1b is given as Remark 3.1 and q denotes the spectral radius of corresponding iteration methods with the above corresponding parameters. The UBOR-type method in Tables 6 and 7 distinguishes the TP-UBOR-type method because of adopting the parameter s different from the above sb : In fact, the value of s is slightly smaller than the one of sb in our experiments. We list the optimal parameters and the corresponding spectral radii in Tables 2 and 4, and also compare the convergence results in Tables 3 and 5 when the three iteration methods with the optimal parameter(s) are employed to solve Problem
Table 1 Values of m; n nonzeros in A and B for Problem 3 Problem 3
m
n
nzðAÞ
nzðBÞ
Stokes 1 Stokes 2 Stokes 3
578 2178 8450
80 288 1088
6178 28,418 122,206
2311 10,453 44,307
669
X.-F. Peng, W. Li / Applied Mathematics and Computation 203 (2008) 660–671 Table 2 Optimal parameters and spectral radius for Problem 1: Q ¼ BT B ðm; nÞ
(50, 40)
(200, 150)
(400, 300)
SOR-like
xb q
1.820069 0.965386
1.953314 0.990364
1.975862 0.995094
GSOR
xb sb q
0.866757 24.071065 0.365024
0.891226 101.675269 0.329808
0.890078 2.016827 0.331544
TP-UBOR-type
sb x2b q
0.866757 0.963992 0.365024
0.891226 0.991235 0.329808
0.890078 0.995587 0.331544
Table 3 The convergence data for Problem 1: Q ¼ BT B ðm; nÞ
(50, 40)
(200, 150)
(400, 300)
SOR-like
IT CPU
314 0.015
1115 0.156
2181 0.453
GSOR
IT CPU
18 0
18 0
18 0
TP-UBOR-type
IT CPU
17 0
16 0.013
16 0.125
Table 4 Optimal parameters and spectral radius for Problem 2 ðm; nÞ Case 1 SOR-like GSOR
TP-UBOR-type
Case 2 SOR-like GSOR
TP-UBOR-type
(128, 64)
(512, 256)
(1152, 576)
xb q xb sb q sb x2b q
1.135249 0.813556 0.797409 2.247650 0.450101 0.797409 0.645226 0.450101
1.152876 0.891923 0.689291 2.983467 0.557413 0.689291 0.768963 0.557413
1.157750 0.923075 0.620632 3.536992 0.615929 0.620632 0.824531 0.615929
xb q xb sb q sb x2b q
1.112822 0.802396 0.799522 2.095872 0.447747 0.799522 0.618525 0.447747
1.133616 0.889414 0.685604 2.843637 0.560710 0.685604 0.758899 0.560710
1.141254 0.922135 0.616097 3.406877 0.619599 0.616097 0.819161 0.619599
Table 5 The convergence data for Problem 2 ðm; nÞ Case 1
SOR-like GSOR TP-UBOR-type
Case 2
SOR-like GSOR TP-UBOR-type
(128, 64)
(512, 256)
(1152, 576)
IT CPU IT CPU IT CPU
63 0.047 22 0.016 21 0.015
120 0.422 31 0.110 30 0.297
173 1.562 38 0.360 38 1.732
IT CPU IT CPU IT CPU
60 0.047 22 0.031 21 0.078
110 1.062 31 0.328 30 2.407
160 1.469 38 0.343 37 4.125
670
X.-F. Peng, W. Li / Applied Mathematics and Computation 203 (2008) 660–671
Table 6 Parameters for Problem 3 Problem 3 Case 1
SOR-like GSOR UBOR-type
Case 2
SOR-like GSOR UBOR-type
Stokes 1
Stokes 2
Stokes 3
xb xb sb s x1b x2b
0.934741 0.074732 21.130564 0.067800 260.620141 0.996463
0.902624 0.034160 42.725887 0.032100 1207.037537 0.999200
0.887112 0.016318 85.960091 0.015650 5180.684141 0.999810
xb xb sb s x1b x2b
0.981697 0.068948 26.506476 0.062500 356.936524 0.997399
0.957138 0.031042 55.719507 0.029200 1738.246344 0.999443
0.940909 0.014625 113.477102 0.014000 7644.632445 0.999871
Table 7 Iteration numbers for Problem 3 Problem 3
Stokes 1
Stokes 2
Stokes 3
Case 1
SOR-like GSOR UBOR-type
7241 379 265
19,261 838 526
29,987 1739 1009
Case 2
SOR GSOR UBOR-type
13,822 455 291
> 30; 000 977 559
> 30; 000 2006 1032
1–2. Evidently, the spectral radii of the TP-UBOR-type and GSOR methods are much smaller than that of the SOR-like method in these cases and the former two methods are much more effective than the latter in the sense of the number of iterations. The numerical results are in agreement with the theoretical results in Theorem 3.1, namely, the TP-UBOR-type method has almost the same convergence rate as the GSOR method when the above optimal parameters are employed. As predicted, the GSOR method takes less CPU time than the TP-UBOR-type with the same convergence rate because of its less computational cost than the latter at each iteration step. However, as showed in Tables 6 and 7, the UBOR-type method can also outperform considerably the GSOR method with optimal parameters in iteration steps, where the parameter s is used instead of sb in the TP-UBOR-type method. Remark 4.1. It is worth mentioning that the relaxation parameters of the UBOR-type method in Tables 6 and 7 are not optimal and only lie in the convergence region of the method. Therefore, the UBOR-type method with the optimal parameters will be considerably effective due to the least iteration steps and in the wake of the significant decrease on the CPU time.
5. Conclusions and remarks In the paper, we have put forward a new unsymmetric block overrelaxation-type (UBOR-type) method for solving the saddle point problem, which is a generalization of the SOR-like method and is completely different from the GSOR method in [22] in algorithmic structure. In theory, we have investigated the convergence and the optimal parameters of the new method under the restriction s ¼ x1 þ x2 x1 x2 . It has been proven that the UBOR-type method has much faster convergence rate than the SOR-like method. Moreover, the new method can achieve the same asymptotic convergence rate as the GSOR method at least. In practice, our method may also be much more effective than the GSOR method in the sense of iteration step as showed in Table 7. In future, we shall pay great attention to determining the optimal parameters of the UBORtype method without the restriction s ¼ x1 þ x2 x1 x2 . In addition, we shall also further optimize our algorithm to reduce the computational cost. Acknowledgement The authors would like to thank the anonymous referee for his/her helpful suggestions which greatly improve this paper. References [1] F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, Springer-Verlag, New York, 1991. [2] S. Wright, Stability of augmented system factorizations in interior-point methods, SIAM J. Matrix Anal. Appl. 18 (1997) 191–222.
X.-F. Peng, W. Li / Applied Mathematics and Computation 203 (2008) 660–671 [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24]
671
B. Fischer, A. Ramage, D.J. Silvester, A.J. Wathen, Minimum residual methods for augmented systems, BIT 38 (1998) 527–543. H.C. Elman, G.H. Golub, Inexact and preconditioned Uzawa algorithms for saddle point problems, SIAM J. Numer. Anal. 31 (1994) 1645–1661. 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. W. Zulehner, Analysis of iterative methods for saddle point problems: a unified approach, Math. Comput. 71 (2002) 479–505. Q.-Y. Hu, J. Zou, An iterative method with variable relaxation parameters for saddle-point problems, SIMA J. Matrix Anal. Appl. 23 (2001) 317–338. G.H. Golub, A.J. Wathen, An iteration for indefinite systems and its application to the Navier–Stokes equations, SIAM J. Sci. Comput. 19 (1998) 530–539. G.H. Golub, X. Wu, J.Y. Yuan, SOR-like methods for augmented systems, BIT 41 (2001) 71–85. 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. 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. M. Benzi, G.H. Golub, A preconditioner for generalized saddle point problems, SIAM J. Matrix Anal. Appl. 26 (2004) 20–41. Z.-Z. Bai, G.H. Golub, L.-Z. Lu, J.-F. Yin, Block triangular and skew-Hermitian splitting methods for positive-definite linear systems, SIAM J. Sci. Comput. 26 (2005) 844–863. 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. C. Keller, N.I.M. Golud, A.J. Wathen, Constraint preconditioning for indefinite linear systems, SIAM J. Matrix Anal. Appl. 21 (2000) 1300–1317. E. DeSturler, J. Liensen, Block-diagonal and constraint preconditioners for nonsymmetric indefinite linear systems. Part: Theory, SIAM J. Sci. Comput. 26 (2005) 1598–1619. Z.-Z. Bai, G.-Q. Li, Restrictively preconditioned conjugate gradient methods for systems of linear equations, IMA J. Numer. Anal. 23 (2003) 561–580. Z.-Z. Bai, Z.-Q. Wang, Restrictive preconditioners for conjugate gradient methods for symmetric positive definite linear systems, J. Comput. Appl. Math. 187 (2006) 202–226. J.-Y. Pan, M.K. Ng, Z.-Z. Bai, New preconditioners for saddle point problems, Appl. Math. Comput. 172 (2006) 762–771. Z. Li, Iterative methods for large sparse saddle point linear systems, Ph.D. Thesis, University of Northeastern, China, 2005 (in Chinese). X.-H. Shao, H.-L. Shen, C.-J. Li, T. Zhang, Generalized AOR method for augmented systems, J. Numer. Meth. Comput. Appl. 27 (2006) 241–248 (in Chinese). Z.-Z. Bai, B.N. Parlett, Z.-Q. Wang, On generalized successive overrelaxation methods for augmented linear systems, Numer. Math. 102 (2005) 1–38. D.M. Young, Iterative Solutions of Large Linear Systems, Academic Press, New York, 1971. R.A. Horn, C.R. Johnson, Matrix Analysis, Cambridge University Press, New York, 1985.
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
On unsymmetric block overrelaxation-type methods for saddle point problems q Xiao-Fei Peng, Wen Li * School of Mathematical Sciences, South China Normal University, Guangzhou 510631, PR China
a r t i c l e
i n f o
a b s t r a c t The unsymmetric block overrelaxation-type (UBOR-type) method is proposed to attack saddle point problems in this paper. The convergence and the optimal parameters for the method are studied when the iteration parameters satisfy some relationship. Theoretical analyses show that the UBOR-type method has faster asymptotic convergence rate than the SOR-like method and its convergence rate can reach the same as that of the GSOR method at least. Numerical experiments support our theoretical results. Moveover, the numerical results further reveal that the new method can be much more effective than the GSOR method in terms of iteration steps. Ó 2008 Elsevier Inc. All rights reserved.
Keywords: The saddle point problem The UBOR-type method Convergence Optimal parameters
1. Introduction Consider the saddle point problem
A
B T
B
0
x y
¼
b q
;
or
eu ~ ¼ ~f ; A
ð1:1Þ
where A 2 Rmm is a symmetric positive definite matrix, B 2 Rmn is a matrix of full column rank with m > n and b 2 Rm , q 2 Rn are two given vectors. The saddle point problem appears frequently in many scientific and engineering computing [1–3]. A number of solution methods for solving the linear system (1.1) can be found in the literature, such as Uzawa-type schemes [4–7], splitting methods [8–14], preconditioned Krylov subspace methods [12,15–19]. Recently, Golub et al. proposed the SOR-like method and considered the optimum choice for the iterative parameter by using some nonsingular preconditioning e [9]. Li and Shao et al. generalized the relaxed splitting method to the matrix Q instead of the null block in coefficient matrix A GAOR method to improve the convergence rate [20,21]. However, they do not achieve the optimal parameters, which affects the efficiency of the method. By using a parameter matrix instead of a single parameter, Bai et al. presented the GSOR method and accelerated the SOR-like method without extra cost per iteration step [22]. Additionally, by using the techniques of vector extrapolation, matrix relaxation and inexact iteration, the GSOR method was extended to the generalized inexact accelerated overrelaxation (GIAOR) method in [22]. To further improve the convergence rate, we present an unsymmetric block overrelaxation-type (UBOR-type) method with three parameters s; x1 ; x2 in this paper. In Section 2, the new splitting method is proposed and its convergence properties are studied. In Section 3, the optimal parameters and corresponding spectral radius are determined when the param-
q
The work was supported by National Natural Science Foundation (Grant No. 10671077) and Guangdong provincial Natural Science Foundation (Grant No. 06025061). * Corresponding author. E-mail addresses: [email protected] (X.-F. Peng), [email protected] (W. Li).
0096-3003/$ - see front matter Ó 2008 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2008.05.014
661
X.-F. Peng, W. Li / Applied Mathematics and Computation 203 (2008) 660–671
eters satisfy the relationship s ¼ x1 þ x2 x1 x2 . Numerical experiments are presented to support our theoretical results in Section 4. Finally, we draw our conclusions in Section 5. 2. The UBOR-type method Let Q 2 Rnn be nonsingular and symmetric, and
e¼ A
A
B
BT
0
¼ D L U;
where
D¼
A 0 ; 0 Q
L¼
0 ; 0
0 BT
U¼
0 B : 0 Q
Then the UBOR-type iteration scheme for solving the augmented linear system (1.1) is given by
xðkþ1Þ
! ¼ Gs;x1 ;x2
yðkþ1Þ
xðkÞ
! þ M 1 s;x1 ;x2
yðkÞ
b ; q
and the iteration matrix Gs;x1 ;x2 of the UBOR-type method can be written as
Gs;x1 ;x2 ¼ M1 s;x1 ;x2 N s;x1 ;x2 ; where
M s;x1 ;x2 ¼
1
s
ðD x1 LÞD1 ðD x2 UÞ ¼
1
A T
s x1 B
x2 B ð1 x2 ÞQ x1 x2 BT A1 B
and
e¼ Ns;x1 ;x2 ¼ M s;x1 ;x2 A ¼
1
s
ð1 sÞA s ðs x1 ÞBT
1
½ð1 sÞD þ ðs x1 ÞL þ ðs x2 ÞU þ x1 x2 LD1 U ! ðx2 sÞB : ð1 x2 ÞQ x1 x2 BT A1 B
Using the block-inversion formula
A11
A12
A21
A22
1 ¼
1 1 1 ½I þ A1 11 A12 ðA22 A21 A11 A12 Þ A21 A11
1 1 A1 11 A12 ðA22 A21 A11 A12 Þ
1 1 ðA22 A21 A1 11 A12 Þ A21 A11
1 ðA22 A21 A1 11 A12 Þ
! ;
we obtain
0 ðsMs;x1 ;x2 Þ1 ¼ @
1 x2 A1 BQ 1 BT A1 A1 x 1x2
x1
1x2
Q 1 BT A1
1xx2 2 A1 BQ 1 1 1x2
Q 1
1 A:
Hence 1 Gs;x1 ;x2 ¼ M1 s;x ;x N s;x1 ;x2 ¼ ðsM s;x1 ;x2 Þ ðsN s;x1 ;x2 Þ 0 1 2 1 1 1 T 1 T 1 x1 Þx2 s 1 1 x2 s sA1 B þ x1 ð1 sÞI ð11 x2 A BQ B x2 A BQ B A B A @ ; ¼ 1 T 1 ð1x1 Þs 1 T x1 s Q B I 1 1x2 x2 Q B A B
where I is the identity matrix, parameters s; x1 ; x2 are real numbers with
ð2:1Þ
s 6¼ 0; x2 6¼ 1.
Algorithm 2.1 (The UBOR-type method) 1. Given initial vectors xð0Þ ; yð0Þ , the relaxation parameters s; x1 ; x2 , and the preconditioning matrix Q, which is an approximate matrix of the Schur complement matrix BT A1 B. 2. For k ¼ 0; 1; . . . ; until convergence, do 1 T 1 1 ðkÞ T ðkÞ sx1 s yðkþ1Þ ¼ yðkÞ þ 1 x2 Q B A ðb By Þ þ 1x2 Q ½ð1 x1 ÞB x q;
xðkþ1Þ ¼ ð1 sÞxðkÞ þ sA1 ðb ByðkÞ Þ x2 A1 Bðyðkþ1Þ yðkÞ Þ: It is worth mentioning that we have computed yðkþ1Þ before computing xðkþ1Þ . Obviously, when x2 ¼ 0, the above algorithm reduces to the GAOR method [20]:
662
X.-F. Peng, W. Li / Applied Mathematics and Computation 203 (2008) 660–671
xðkþ1Þ ¼ ð1 sÞxðkÞ þ sA1 ðb ByðkÞ Þ; yðkþ1Þ ¼ yðkÞ Q 1 fsq þ BT ½ðx1 sÞxðkÞ x1 xðkþ1Þ g: 1s When three parameters satisfy the relationship s ¼ x1 þ x2 x1 x2 , i.e., x1 ¼ 1 1 x2 , the UBOR-type method reduces to the following two-parameter UBOR-type (TP-UBOR-type) method:
sðsx2 Þ 1 T 1 yðkþ1Þ ¼ yðkÞ þ ð1 Q B A ðb ByðkÞ Þ þ 1sx2 Q 1 x Þ2
2
1s 1x2
BT xðkÞ q ;
xðkþ1Þ ¼ ð1 sÞxðkÞ þ sA1 ðb ByðkÞ Þ x2 A1 Bðyðkþ1Þ yðkÞ Þ: e namely, Remark 2.1. The TP-UBOR-type method is a fixed point iteration method based on the splitting of A, e ¼ Ms;x N s;x , where A 2 2
M s;x2 ¼
!
x2 B
A
1
ðsx2 Þ T 1 sx2 T ð1 x2 ÞQ x21 s 1 x2 B x2 B A B
and
Ns;x2 ¼
ð1 sÞA
1
!
ðx2 sÞB
ðsx2 Þ T 1 2 ð1sÞ T s x1 ð1 x2 ÞQ x21 x2 B x2 B A B
:
Specially, if x2 ¼ 0, the TP-UBOR-type algorithm reduces to the SOR-like method [9]:
xðkþ1Þ ¼ ð1 sÞxðkÞ þ sA1 ðb ByðkÞ Þ; yðkþ1Þ ¼ yðkÞ þ sQ 1 ðBT xðkþ1Þ qÞ: Remark 2.2. It is noted that the TP-UBOR-type method and GSOR method in [22] have different algorithmic structures based on different splittings of the coefficient matrix, but they can both become the SOR-like method under the different restrictions on the iteration parameters. Generally speaking, the computational cost of the former can be more than that of the latter at each step. To improve the running speed, we consider the following practical Algorithm. Algorithm 2.2 (The practical UBOR-type method) 1. Given initial vectors xð0Þ ; yð0Þ , the relaxation parameters s; x1 ; x2 , and the preconditioning matrix Q, which is an approximate matrix of the Schur complement matrix BT A1 B. 2. Let t1 ¼ 1sx2 ; t 0 ¼ t1 x1 ; t 2 ¼ 1 s; t3 ¼ 1 x1 . 3. Let T 1 ¼ Q 1 BT ; T 2 ¼ A1 b; T 3 ¼ A1 B; T 4 ¼ Q 1 q. 4. For k ¼ 0; 1; . . . ; until convergence, do
T ¼ T 2 T 3 yðkÞ ; yðkþ1Þ ¼ yðkÞ þ t0 T 1 T þ t 1 ½t 3 T 1 xðkÞ T 4 ; xðkþ1Þ ¼ t2 xðkÞ þ ðs x2 ÞT þ x2 ðT 2 T 3 yðkþ1Þ Þ: Next, we shall study the convergence of the UBOR-type method for soving the linear system (1.1). By rðÞ and qðÞ we denote the spectrum and spectral radius of a matrix, respectively. Let k 2 rðGs;x1 ;x2 Þ and l 2 rðQ 1 BT A1 BÞ; lmin ¼ minflg and q ¼ qðQ 1 BT A1 BÞ. Suppose that ðxT ; yT ÞT is the eigenvector corresponding to k. Then
Gs;x1 ;x2
x y
x ¼k ; y
i.e.,
ð1 sÞA
ðs x2 ÞB
ðs x1 ÞBT
ð1 x2 ÞQ x1 x2 BT A1 B
! x y
¼k
A x1 BT
x2 B ð1 x2 ÞQ x1 x2 BT A1 B
x ; y
ð2:2Þ
from which one may deduce that
ð1 s kÞx ¼ ðkx2 þ s x2 ÞA1 By
ð2:3Þ
ð1 kÞ½x1 x2 BT A1 By ð1 x2 ÞQy ¼ ðkx1 þ s x1 ÞBT x:
ð2:4Þ
and
X.-F. Peng, W. Li / Applied Mathematics and Computation 203 (2008) 660–671
663
Lemma 2.1 [23]. Both roots of the real quadratic equation x2 bx þ c ¼ 0 are less than one in modulus if and only if jcj < 1 and jbj < 1 þ c: Lemma 2.2. Assume that 1 s 6¼ 0. Then (1) when x1 ¼ 1; k ¼ 1 s is an eigenvalue with m multiplicities at least; (2) when x1 6¼ 1 and m > n, k ¼ 1 s is an eigenvalue with m n multiplicities at least. Proof (1) When x1 ¼ 1, by (2.1) we obtain
0 Gs;x1 ;x2 ¼ @
ð1 sÞI
1 T 1 x2 s 1 sA1 B þ 1 x2 A BQ B A B
O
I 1sx2 Q 1 BT A1 B
1 A;
which implies the conclusion (1). (2) From (2.3) and (2.4), an explicit calculation reveals that k ¼ 1 s is an eigenvalue of Gs;x1 ;x2 with m n multiplicities at least and the corresponding eigenvector is ðxT ; 0ÞT , where ð1 x1 ÞBT x ¼ 0: In terms of the assumptions on x1 and B, it holds that x 2 NðBT Þ and x 6¼ 0. h Theorem 2.1. Let A 2 Rmm and Q 2 Rnn be symmetric positive definite, and B 2 Rmn be of full column rank. Then the UBOR-type method is convergent if and only if
0 < s < 2; x2 < 1
qs2
; 4 s x2 1 s 2x2 2 s < x1 < : þ sq 1 x2 q 2ð1 x2 Þ 1
1
Proof. Since Q 1 BT A1 B is similar to Q 2 BT A1 BQ 2 ; they have same eigenvalues. Moreover, the latter is congruent to the positive definite matrix BT A1 B: It follows from Sylvester’s law of inertia that l > 0 [24]. By Lemma 2.2, k ¼ 1 s 6¼ 0 is an eigenvalue of Gs;x1 ;x2 : If k 6¼ 1 s, it follows from (2.3) and (2.4) that
k2 ðsl slx1 þ 2 s
sl 1 x2
Þk þ sl slx1 þ 1 s
slð1 sÞ ¼ 0; 1 x2
ð2:5Þ
or equivalently,
k2 ðc þ 1
s2 l Þk þ c ¼ 0; 1 x2
where
c ¼ sl slx1 þ 1 s
slð1 sÞ : 1 x2
By Lemma 2.1, jkj < 1 if and only if
j1 sj < 1; jcj < 1;
2
s l jc þ 1 1 x2 j < 1 þ c:
It follows that
0 < s < 2;
x2 < 1;
and
s2 l 2ð1 x2 Þ
< c þ 1 < 2;
s2 l 2ð1 x2 Þ
< 2:
In terms of (2.2), we have
s x2 1 s x2 1 6 < x1 ; 1 x2 l 1 x2 q s 2x2 2 s s 2x2 2 s x1 < 6 : þ þ sq sl 2ð1 x2 Þ 2ð1 x2 Þ s2 l 6 s2 q < 4ð1 x2 Þ;
ð2:6Þ
664
X.-F. Peng, W. Li / Applied Mathematics and Computation 203 (2008) 660–671
Therefore, we obtain the following necessary and sufficient conditions such that the UBOR-type method is convergent
0 < s < 2; x2 < 1
qs2
; 4 s x2 1 s 2x2 2 s < x1 < : þ sq 1 x2 q 2ð1 x2 Þ
Let x2 ¼ 0; then the following result follows immediately from Theorem 2.1. Corollary 2.1 [20]. Under the same assumptions as in Theorem 2.1, the GAOR method is convergent if and only if
1 s 2s s < x1 < þ : sq q 2
2 0 < s < minf2; pffiffiffiffig;
q
Remark 2.3. Actually, the above GAOR method is a special GIAOR method in [22] since the latter may reduce to the former when we take P ¼ A, D ¼ cI and X ¼ xI. If
s ¼ x1 þ x2 x1 x2 , by Theorem 2.1 then we have:
Corollary 2.2. Under the same assumptions as in Theorem 2.1, the TP-UBOR-type method is convergent if and only if
0 < s < 2;
x2 < 1
q s2 2ð2 sÞ
:
Corollary 2.3 [9]. Under the same assumptions as in Theorem 2.1, the SOR-like method is convergent if and only if
4 0 < s < pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : 4q þ 1 þ 1 3. The optimal iteration parameters By the proof of Theorem 2.1, k ¼ 1 s or
k1;2
1 ¼ 2
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi)
(
sl l 2 4l ; sl slx1 þ 2 s s ðlx1 þ 1 l þ Þ 1 x2 1 x2 1 x2
ð3:1Þ
where k1 and k2 correspond to taking the ‘+’ and ‘’ signs in (3.1), respectively. Then, we have
qðGs;x1 ;x2 Þ ¼ max fj1 sj; jk1 j; jk2 jg: l According to Theorem 2.1, the parameters of the convergent UBOR-type method satisfy the following conditions:
0 < s < 2; x2 < 1
qs2
; 4 s x2 1 s 2x2 2 s < x1 < : þ sq 1 x2 q 2ð1 x2 Þ
To achieve the absolute values of the nonzero eigenvalues of the iteration matrix Gs;x1 ;x2 , we consider the following three cases: (1)
sl slx1 þ 2 s
(2)
sl slx1 þ 2 s
(3) ðlx1 þ 1 l þ
sl 1 x2
sl 1 x2
l 1 x2
Þ2
P 0;
ðlx1 þ 1 l þ
6 0;
ðlx1 þ 1 l þ
l 1 x2
l 1 x2
Þ2
4l P 0; 1 x2
Þ2
4l P 0; 1 x2
4l 6 0: 1 x2
Accordingly, we define the following functions, whose variables s; x1 ; x2 conform to the restrictions of the above three cases, respectively.
X.-F. Peng, W. Li / Applied Mathematics and Computation 203 (2008) 660–671
1 f1 ðs; x1 ; x2 ; lÞ ¼ 2 for
665
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi)
(
sl l 2 4l sl slx1 þ 2 s þ s ðlx1 þ 1 l þ Þ ; 1 x2 1 x2 1 x2
2 1 x2 1 2 x2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 6 x1 6 ð 1Þ ; l 1 x l s 1 x2 lð1 x2 Þ 2
s6
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 x2 pffiffiffiffi ;
l
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi) ( 1 sl l 2 4l f2 ðs; x1 ; x2 ; lÞ ¼ slx1 sl þ s 2 þ þ s ðlx1 þ 1 l þ Þ ; 2 1 x2 1 x2 1 x2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 x2 1 2 x2 1 ; s 6 pffiffiffiffi or for x1 P l s 1 x2 l pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 x2 2 1 x2 x1 P pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; s P pffiffiffiffi ; l lð1 x2 Þ l 1 x2 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi slð1 sÞ ; f3 ðs; x1 ; x2 ; lÞ ¼ sl slx1 þ 1 s 1 x2
2 1 x2 for x1 6 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : lð1 x2 Þ l 1 x2
Obviously,
qðGs;x1 ;x2 Þ ¼ j1 sj; or
8 2 ffi 1 x2 ; max f3 ðs; x1 ; x2 ; lÞ; for x1 6 pffiffiffiffiffiffiffiffiffiffiffiffiffi > > lð1x2 Þ l 1x2 l > > > > > 2 > ffi 1 x2 6 x1 6 l1 ð2s 1Þ 1xx2 ; max f1 ðs; x1 ; x2 ; lÞ; for pffiffiffiffiffiffiffiffiffiffiffiffiffi > > 2 lð1x2 Þ l 1x2 l > > > pffiffiffiffiffiffiffiffiffi > < 1x2 s 6 pffiffilffi ; qðGs;x1 ;x2 Þ ¼ > pffiffiffiffiffiffiffiffiffi > > 1x > > max f2 ðs; x1 ; x2 ; lÞ; for x1 P l1 ð2s 1Þ 1xx2 2 ; s 6 pffiffilffi 2 > > l > > > pffiffiffiffiffiffiffiffiffi > > > 1x2 2 > ffi l1 1xx2 ; s P p ffiffilffi : : or x1 P pffiffiffiffiffiffiffiffiffiffiffiffiffi 2
ð3:2Þ
lð1x2 Þ
Specially, we shall determine the optimal iteration parameters of the TP-UBOR-type method. Naturally, the notations
qðGs;x2 Þ and fj ðs; x2 ; lÞ in the TP-UBOR-type method can be used instead of qðGs;x1 ;x2 Þ and fj ðs; x1 ; x2 ; lÞ (j = 1, 2, 3) in the UBOR-type method, respectively. Theorem 3.1. Let A 2 Rmm and Q 2 Rnn be symmetric positive definite, B 2 Rmn be of full column rank. Then the spectral radius of the iteration matrix Gs;x2 of the TP-UBOR-type method is given as follows: (1) for 0 < s < 1;
8 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi s2 lmin slmin 2 4lmin 1 > > > 2 ½2 s 1x2 þ s ð1 þ 1x2 Þ 1x2 ; for > < pffiffiffiffiffiffiffiffiffiffiffiffi qðGs;x2 Þ ¼ for 1 s; > > qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > > : 1 ½s 2 þ s2 q þ s ð1 þ sq Þ2 4q ; for 1x2 1x2 1x2 2
x2 < x2ð1Þ ðsÞ; ð2Þ xð1Þ 2 ðsÞ 6 x2 6 x2 ðsÞ;
ð3:3Þ
2
s q xð2Þ 2 ðsÞ < x2 < 1 2ð2sÞ ;
where
xð1Þ 2 ðsÞ ¼ 1
s2 lmin
pffiffiffiffiffiffiffiffiffiffiffiffi ; ð1 1 sÞ2
x2ð2Þ ðsÞ ¼ 1
s2 q
pffiffiffiffiffiffiffiffiffiffiffiffi : ð1 þ 1 sÞ2
(2) for 1 6 s < 2; we have
8 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 slmin 2 4lmin > < 12 ½2 s s1lxmin þ s ð1 þ 1 x2 Þ 1x2 ; for 2 qðGs;x2 Þ ¼ q ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi > : 1 ½s 2 þ s2 q þ s ð1 þ sq Þ2 4q ; for 1x2 1x2 1x2 2 where
xð3Þ 2 ðsÞ ¼ 1
s2 ðq þ lmin Þ : 2ð2 sÞ
x2 < x2ð3Þ ðsÞ; 2
s q xð3Þ 2 ðsÞ 6 x2 < 1 2ð2sÞ ;
ð3:4Þ
666
X.-F. Peng, W. Li / Applied Mathematics and Computation 203 (2008) 660–671
Furthermore, the optimal parameters are given by
pffiffiffiffiffiffiffiffiffiffiffiffiffi
4 qlmin 4qlmin sb ¼ pffiffiffiffi pffiffiffiffiffiffiffiffiffi ffi 2 ; x2b ¼ 1 pffiffiffiffi pffiffiffiffiffiffiffiffiffiffi 2 ; ð q þ lmin Þ ð q þ lmin Þ and the corresponding optimal spectral radius of the iteration matrix Gs;x2 is
pffiffiffiffi pffiffiffiffiffiffiffiffiffiffi q lmin min qðGs;x2 Þ ¼ qðGsb ;x2b Þ ¼ pffiffiffiffi pffiffiffiffiffiffiffiffiffi ffi: s;x 2 q þ lmin
Proof. By Corollary 2.2, it need only to consider that parameters s and x2 satisfy
0 < s < 2;
x2 < 1
q s2 2ð2 sÞ
:
According to the expressions of the three functions fj ðs; x1 ; x2 ; lÞ ðj ¼ 1; 2; 3Þ we derive the following expressions of the functions fj ðs; x2 ; lÞ ðj ¼ 1; 2; 3Þ in the TP-UBOR-type method as follows:
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffio n s2 l sl 2 4l f1 ðs; x2 ; lÞ ¼ 12 2 s 1 x2 þ s ð1 þ 1x2 Þ 1x2 ; 2
l ; for 0 < s < 1; x2 < 1 ð1ps ffiffiffiffiffiffi 1sÞ2 2
s l or 1 6 s < 2; x2 < 1 2 s; qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffio n 2l s f2 ðs; x2 ; lÞ ¼ 12 s 2 þ 1x2 þ s ð1 þ 1slx2 Þ2 14lx2 ; 2
2
l s q < x2 < 1 2ð2 for 0 < s < 1; 1 ð1þps ffiffiffiffiffiffi sÞ ; 1sÞ2 2
2
s l s q or 1 6 s < 2; 1 2 s 6 x2 < 1 2ð2sÞ ; pffiffiffiffiffiffiffiffiffiffiffiffi f3 ðs; x2 ; lÞ ¼ 1 s; 2
2
l l 6 x2 6 1 ð1þps ffiffiffiffiffiffi : for 0 < s < 1; 1 ð1ps ffiffiffiffiffiffi 1sÞ2 1sÞ2
Under the above restrictions on s; x2 , it holds that f3 ðs; x2 ; lÞ ¼ fj ðs; x2 ; lÞ ðj ¼ 1; 2; 3Þ to achieve qðGs;x2 Þ. A simple computation reveals that
8 > > <
pffiffiffiffiffiffiffiffiffiffiffiffi 1 s P j1 sj: Hence, we just consider the functions
9 > > = of1 ðs; x2 ; lÞ s ¼ s þ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi < 0; > 2 ol 2ð1 x2 Þ > > sl 4l > : ; 1 þ 1x2 1x2 8 9 > > > > s2 l < = s 2 þ 1x2 of2 ðs; x2 ; lÞ s ¼ s þ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > 0; 2 > ol 2ð1 x2 Þ > > > : 1 þ sl 4l ; s2 l 2 s 1 x2
1x2
ð3:5Þ
1x2
of3 ðs; x2 ; lÞ ¼ 0; ol and
8 > > <
9 > > =
s2 l 2 s 1 x2
of1 ðs; x2 ; lÞ sl ffi < 0; ¼ s þ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 > o x2 2ð1 x2 Þ2 > > > : 1 þ sl 4l ; 8 > > <
1x2
1x2
9 > > =
2
s 2 þ 1s xl2
of2 ðs; x2 ; lÞ sl ffi > 0; ¼ s þ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 > o x2 2ð1 x2 Þ2 > > > : 1 þ sl 4l ; 1x2
1x2
of3 ðs; x2 ; lÞ ¼ 0: o x2 When 0 < s < 1; or 1 6 s < 2; then (3.3) or (3.4) follows immediately from (3.2) and (3.5). Now, we consider the optimal relaxation parameters in the following two cases. (1) For the case 0 < s < 1; it follows from (3.3) and (3.6) that
ð3:6Þ
667
X.-F. Peng, W. Li / Applied Mathematics and Computation 203 (2008) 660–671
8 ð1Þ ð1Þ > min f1 ðs; x2 ; lmin Þ ¼ f1 ðs; x2 ðsÞ; lmin Þ; x2 6 x2 ðsÞ; > > > x2 < pffiffiffiffiffiffiffiffiffiffiffiffi ð2Þ min qðGs;x2 Þ ¼ xð1Þ 1 s; 2 ðsÞ 6 x2 6 x2 ðsÞ; x2 > > > ð2Þ ð2Þ s2 q > : min f2 ðs; x2 ; qÞ ¼ f2 ðs; x2 ðsÞ; qÞ; x2 ðsÞ 6 x2 < 1 2ð2 sÞ :
ð3:7Þ
x2
For any fixed s, it is easy to verify that ð1Þ
ð2Þ
f1 ðs; x2 ðsÞ; lmin Þ ¼ f2 ðs; x2 ðsÞ; qÞ ¼
pffiffiffiffiffiffiffiffiffiffiffiffi 1 s: ð1Þ
Evidently, minx2 qðGs;x2 Þ attains the minimum when the parameter s attains the maximum. Here, it noted that x2 ðsÞ and ðiÞ xð2Þ 2 ðsÞ are the functions of s. Next, we investigate further the monotonicity of x2 ðsÞ; i ¼ 1; 2. Because
pffiffiffiffiffiffiffiffiffiffiffiffi ð1Þ ox2 ðsÞ lmin ð1 þ 1 sÞ pffiffiffiffiffiffiffiffiffiffiffiffi > 0; ¼ os 1s
pffiffiffiffiffiffiffiffiffiffiffiffi ð2Þ ox2 ðsÞ qð1 1 sÞ pffiffiffiffiffiffiffiffiffiffiffiffi < 0; ¼ os 1s
ð2Þ xð1Þ 2 ðsÞ is increasing and x2 ðsÞ is decreasing when s is increasing. To ensure that s both attains the maximum and satisfies ð1Þ
ð2Þ
ð1Þ
ð2Þ
the inequality x2 ðsÞ 6 x2 ðsÞ, the s such that qðGs;x2 Þ attains the minimum must satisfy the relation x2 ðsÞ ¼ x2 ðsÞ; namely,
1
s2 lmin
s2 q
pffiffiffiffiffiffiffiffiffiffiffiffi ¼ 1 pffiffiffiffiffiffiffiffiffiffiffiffi : ð1 1 sÞ2 ð1 þ 1 sÞ2
It then follows that ð2Þ xð1Þ qðGs;x2 Þ 2 ðsÞ ¼ x2 ðsÞ ¼ arg min x 2
for the fixed s. An simple calculation shows that the optimal relaxation parameters are
pffiffiffiffiffiffiffiffiffiffiffiffiffi
4 qlmin 4qlmin sbð1Þ ¼ pffiffiffiffi pffiffiffiffiffiffiffiffiffi ffi ; x2bð1Þ ¼ 1 pffiffiffiffi pffiffiffiffiffiffiffiffiffiffi 2 ; ð q þ lmin Þ ð q þ lmin Þ2 and the corresponding optimal spectral radius is
pffiffiffiffi
pffiffiffiffiffiffiffiffiffiffi
q lmin qðGðsbð1Þ ; x2bð1Þ ÞÞ ¼ pffiffiffiffi pffiffiffiffiffiffiffiffiffi ffi: q þ lmin (2) For the case 1 6 s < 2; it follows from (3.4) and (3.6) that ð3Þ
ð3Þ
min qðGs;x2 Þ ¼ f1 ðs; x2 ðsÞ; lmin Þ ¼ f2 ðs; x2 ðsÞ; qÞ x2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 ¼ ðq lmin Þð2 sÞ þ s2 ðq þ lmin Þ2 4qlmin ð2 sÞ2 : 2ðq þ lmin Þ Let yðsÞ ¼ 2ðq þ lmin Þ minx2 qðGðs; x2 ÞÞ: Because
oy ðq lmin Þ2 s þ 8qlmin P 0; ¼ ðq lmin Þ þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi os s2 ðq þ lmin Þ2 4qlmin ð2 sÞ2 the optimal relaxation parameters such that qðGðs; x2 ÞÞ attains the minimum are
sbð2Þ ¼ 1; x2bð2Þ ¼ xð3Þ 2 ¼ 1
q þ lmin 2
;
and the corresponding optimal spectral radius is
qðGsbð2Þ ;x2bð2Þ Þ ¼
q lmin : q þ lmin
Comparing the above two cases, the desired results follow. h Remark 3.1. It is noted that the sb of the optimal TP-UBOR-type method is just equal to the xb of the optimal GSOR method in [22] and their corresponding methods have the same asymptotic convergence factor. In other words, the UBOR-type method is equally effective to the GSOR method at least in theory since the TP-UBOR-type method is a special UBOR-type method.
668
X.-F. Peng, W. Li / Applied Mathematics and Computation 203 (2008) 660–671
Remark 3.2. Generally speaking, the optimal parameters sb ; x2b of the TP-UBOR-type method and the corresponding x1b are not optimal in the UBOR-type method, where the parameter x1b is determined by the relationship sb ¼ x1b þ x2b x1b x2b . By Theorem 3.1,
x1b ¼ 1
pffiffiffiffi
pffiffiffiffiffiffiffiffiffiffi
q lmin 2 : 4qlmin
It is very difficult to determine the optimal relaxation parameters of the UBOR-type method without the above given relationship on three parameters. However, as we shall show in next section on numerical experiments, the faster convergence rate of the UBOR-type method can be achieved for other cases different from those in Theorem 3.1.
4. Numerical experiments In this section, we apply the UBOR-type method to solve those problems taken from [7] and [22], and a set of Stokes saddle point problems generated using the IFISS software written by Howard Elman et al. Problem 1. Consider the linear system (1.1) with the coefficient sub-matrices A; B deriving from [7],
8 > < i þ 1; i ¼ j; A ¼ 1; ji jj ¼ 1; > : 0; other;
B¼
j;
i ¼ j þ m n;
0; other:
The right vectors are given by
b ¼ ð1; 1; . . . ; 1ÞT 2 Rm ;
q ¼ ð1; 1; . . . ; 1ÞT 2 Rn :
Problem 2. Consider the linear system with the coefficient sub-matrices A; B in [22],
A¼
IT þT I
0
0
IT þT I
;
B¼
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. The right vectors are defined as
b ¼ ð1; 1; . . . ; 1ÞT 2 Rm ;
q ¼ ð0; 0; . . . ; 0ÞT 2 Rn ; m ¼ 2p2 ; n ¼ p2 :
Problem 3. We generate the regularised-lid driven cavity problems with IFISS software using the 16 16; 32 32; 64 64 grid, Q 2 Q 1 discretisation. Since the matrix B generated by the package is actually singular, we drop the first two rows for the Stokes 1 problem and the first rows of B for the Stokes 2–3 problems to obtain the full rank matrices. The dimensions of the problems and number of nonzeros in the associated matrices A and B are given in Table 1. All experiments are performed in MATLAB with the machine precision 1016 and all numerical results are listed in Tables 2–7. To save the storage space, we just list the optimal parameters and spectral radius with six number of decimal place in 6 kr k k following Tables. We use a zero initial guess and the stopping criterion kr , where r k ; r 0 denote the kth residual vector 0 k 6 10 and the initial one, respectively. IT and CPU represent the number of iteration steps and the elasped CPU time in seconds, respectively. The preconditioning matrix Q is chosen as BT B for Problem 1 and is given as follows for Problems 2 and 3,
Case 1 : Q ¼ diagðBT A1 BÞ;
Case 2 : Q ¼ tridiagðBT A1 BÞ:
xb ; sb and x2b indicate the corresponding optimal parameters in the SOR-like, GSOR and UBOR-type methods. x1b is given as Remark 3.1 and q denotes the spectral radius of corresponding iteration methods with the above corresponding parameters. The UBOR-type method in Tables 6 and 7 distinguishes the TP-UBOR-type method because of adopting the parameter s different from the above sb : In fact, the value of s is slightly smaller than the one of sb in our experiments. We list the optimal parameters and the corresponding spectral radii in Tables 2 and 4, and also compare the convergence results in Tables 3 and 5 when the three iteration methods with the optimal parameter(s) are employed to solve Problem
Table 1 Values of m; n nonzeros in A and B for Problem 3 Problem 3
m
n
nzðAÞ
nzðBÞ
Stokes 1 Stokes 2 Stokes 3
578 2178 8450
80 288 1088
6178 28,418 122,206
2311 10,453 44,307
669
X.-F. Peng, W. Li / Applied Mathematics and Computation 203 (2008) 660–671 Table 2 Optimal parameters and spectral radius for Problem 1: Q ¼ BT B ðm; nÞ
(50, 40)
(200, 150)
(400, 300)
SOR-like
xb q
1.820069 0.965386
1.953314 0.990364
1.975862 0.995094
GSOR
xb sb q
0.866757 24.071065 0.365024
0.891226 101.675269 0.329808
0.890078 2.016827 0.331544
TP-UBOR-type
sb x2b q
0.866757 0.963992 0.365024
0.891226 0.991235 0.329808
0.890078 0.995587 0.331544
Table 3 The convergence data for Problem 1: Q ¼ BT B ðm; nÞ
(50, 40)
(200, 150)
(400, 300)
SOR-like
IT CPU
314 0.015
1115 0.156
2181 0.453
GSOR
IT CPU
18 0
18 0
18 0
TP-UBOR-type
IT CPU
17 0
16 0.013
16 0.125
Table 4 Optimal parameters and spectral radius for Problem 2 ðm; nÞ Case 1 SOR-like GSOR
TP-UBOR-type
Case 2 SOR-like GSOR
TP-UBOR-type
(128, 64)
(512, 256)
(1152, 576)
xb q xb sb q sb x2b q
1.135249 0.813556 0.797409 2.247650 0.450101 0.797409 0.645226 0.450101
1.152876 0.891923 0.689291 2.983467 0.557413 0.689291 0.768963 0.557413
1.157750 0.923075 0.620632 3.536992 0.615929 0.620632 0.824531 0.615929
xb q xb sb q sb x2b q
1.112822 0.802396 0.799522 2.095872 0.447747 0.799522 0.618525 0.447747
1.133616 0.889414 0.685604 2.843637 0.560710 0.685604 0.758899 0.560710
1.141254 0.922135 0.616097 3.406877 0.619599 0.616097 0.819161 0.619599
Table 5 The convergence data for Problem 2 ðm; nÞ Case 1
SOR-like GSOR TP-UBOR-type
Case 2
SOR-like GSOR TP-UBOR-type
(128, 64)
(512, 256)
(1152, 576)
IT CPU IT CPU IT CPU
63 0.047 22 0.016 21 0.015
120 0.422 31 0.110 30 0.297
173 1.562 38 0.360 38 1.732
IT CPU IT CPU IT CPU
60 0.047 22 0.031 21 0.078
110 1.062 31 0.328 30 2.407
160 1.469 38 0.343 37 4.125
670
X.-F. Peng, W. Li / Applied Mathematics and Computation 203 (2008) 660–671
Table 6 Parameters for Problem 3 Problem 3 Case 1
SOR-like GSOR UBOR-type
Case 2
SOR-like GSOR UBOR-type
Stokes 1
Stokes 2
Stokes 3
xb xb sb s x1b x2b
0.934741 0.074732 21.130564 0.067800 260.620141 0.996463
0.902624 0.034160 42.725887 0.032100 1207.037537 0.999200
0.887112 0.016318 85.960091 0.015650 5180.684141 0.999810
xb xb sb s x1b x2b
0.981697 0.068948 26.506476 0.062500 356.936524 0.997399
0.957138 0.031042 55.719507 0.029200 1738.246344 0.999443
0.940909 0.014625 113.477102 0.014000 7644.632445 0.999871
Table 7 Iteration numbers for Problem 3 Problem 3
Stokes 1
Stokes 2
Stokes 3
Case 1
SOR-like GSOR UBOR-type
7241 379 265
19,261 838 526
29,987 1739 1009
Case 2
SOR GSOR UBOR-type
13,822 455 291
> 30; 000 977 559
> 30; 000 2006 1032
1–2. Evidently, the spectral radii of the TP-UBOR-type and GSOR methods are much smaller than that of the SOR-like method in these cases and the former two methods are much more effective than the latter in the sense of the number of iterations. The numerical results are in agreement with the theoretical results in Theorem 3.1, namely, the TP-UBOR-type method has almost the same convergence rate as the GSOR method when the above optimal parameters are employed. As predicted, the GSOR method takes less CPU time than the TP-UBOR-type with the same convergence rate because of its less computational cost than the latter at each iteration step. However, as showed in Tables 6 and 7, the UBOR-type method can also outperform considerably the GSOR method with optimal parameters in iteration steps, where the parameter s is used instead of sb in the TP-UBOR-type method. Remark 4.1. It is worth mentioning that the relaxation parameters of the UBOR-type method in Tables 6 and 7 are not optimal and only lie in the convergence region of the method. Therefore, the UBOR-type method with the optimal parameters will be considerably effective due to the least iteration steps and in the wake of the significant decrease on the CPU time.
5. Conclusions and remarks In the paper, we have put forward a new unsymmetric block overrelaxation-type (UBOR-type) method for solving the saddle point problem, which is a generalization of the SOR-like method and is completely different from the GSOR method in [22] in algorithmic structure. In theory, we have investigated the convergence and the optimal parameters of the new method under the restriction s ¼ x1 þ x2 x1 x2 . It has been proven that the UBOR-type method has much faster convergence rate than the SOR-like method. Moreover, the new method can achieve the same asymptotic convergence rate as the GSOR method at least. In practice, our method may also be much more effective than the GSOR method in the sense of iteration step as showed in Table 7. In future, we shall pay great attention to determining the optimal parameters of the UBORtype method without the restriction s ¼ x1 þ x2 x1 x2 . In addition, we shall also further optimize our algorithm to reduce the computational cost. Acknowledgement The authors would like to thank the anonymous referee for his/her helpful suggestions which greatly improve this paper. References [1] F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, Springer-Verlag, New York, 1991. [2] S. Wright, Stability of augmented system factorizations in interior-point methods, SIAM J. Matrix Anal. Appl. 18 (1997) 191–222.
X.-F. Peng, W. Li / Applied Mathematics and Computation 203 (2008) 660–671 [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24]
671
B. Fischer, A. Ramage, D.J. Silvester, A.J. Wathen, Minimum residual methods for augmented systems, BIT 38 (1998) 527–543. H.C. Elman, G.H. Golub, Inexact and preconditioned Uzawa algorithms for saddle point problems, SIAM J. Numer. Anal. 31 (1994) 1645–1661. 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. W. Zulehner, Analysis of iterative methods for saddle point problems: a unified approach, Math. Comput. 71 (2002) 479–505. Q.-Y. Hu, J. Zou, An iterative method with variable relaxation parameters for saddle-point problems, SIMA J. Matrix Anal. Appl. 23 (2001) 317–338. G.H. Golub, A.J. Wathen, An iteration for indefinite systems and its application to the Navier–Stokes equations, SIAM J. Sci. Comput. 19 (1998) 530–539. G.H. Golub, X. Wu, J.Y. Yuan, SOR-like methods for augmented systems, BIT 41 (2001) 71–85. 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. 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. M. Benzi, G.H. Golub, A preconditioner for generalized saddle point problems, SIAM J. Matrix Anal. Appl. 26 (2004) 20–41. Z.-Z. Bai, G.H. Golub, L.-Z. Lu, J.-F. Yin, Block triangular and skew-Hermitian splitting methods for positive-definite linear systems, SIAM J. Sci. Comput. 26 (2005) 844–863. 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. C. Keller, N.I.M. Golud, A.J. Wathen, Constraint preconditioning for indefinite linear systems, SIAM J. Matrix Anal. Appl. 21 (2000) 1300–1317. E. DeSturler, J. Liensen, Block-diagonal and constraint preconditioners for nonsymmetric indefinite linear systems. Part: Theory, SIAM J. Sci. Comput. 26 (2005) 1598–1619. Z.-Z. Bai, G.-Q. Li, Restrictively preconditioned conjugate gradient methods for systems of linear equations, IMA J. Numer. Anal. 23 (2003) 561–580. Z.-Z. Bai, Z.-Q. Wang, Restrictive preconditioners for conjugate gradient methods for symmetric positive definite linear systems, J. Comput. Appl. Math. 187 (2006) 202–226. J.-Y. Pan, M.K. Ng, Z.-Z. Bai, New preconditioners for saddle point problems, Appl. Math. Comput. 172 (2006) 762–771. Z. Li, Iterative methods for large sparse saddle point linear systems, Ph.D. Thesis, University of Northeastern, China, 2005 (in Chinese). X.-H. Shao, H.-L. Shen, C.-J. Li, T. Zhang, Generalized AOR method for augmented systems, J. Numer. Meth. Comput. Appl. 27 (2006) 241–248 (in Chinese). Z.-Z. Bai, B.N. Parlett, Z.-Q. Wang, On generalized successive overrelaxation methods for augmented linear systems, Numer. Math. 102 (2005) 1–38. D.M. Young, Iterative Solutions of Large Linear Systems, Academic Press, New York, 1971. R.A. Horn, C.R. Johnson, Matrix Analysis, Cambridge University Press, New York, 1985.