- Email: [email protected]
Semi-convergence analysis of Uzawa methods for singular saddle point problems
Journal of Computational and Applied Mathematics 255 (2014) 334–345
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
Semi-convergence analysis of Uzawa methods for singular saddle point problems Naimin Zhang a , Tzon-Tzer Lu b , Yimin Wei c,d,∗ a
School of Mathematics and Information Science, Wenzhou University, Wenzhou, 325035, PR China
b
Department of Applied Mathematics, National Sun Yat-sen University, Kaohsiung, 80424, Taiwan
c
School of Mathematical Sciences, Fudan University, Shanghai, 200433, PR China
d
Shanghai Key Laboratory of Contemporary Applied Mathematics, PR China
article
abstract
info
Article history: Received 19 January 2012 Received in revised form 28 November 2012 MSC: 15A09 65F10 Keywords: Singular linear systems Saddle point problems Uzawa method Semi-convergence
Recently, Zheng, Bai and Yang studied the parameterized Uzawa method for solving singular saddle point problems (B. Zheng, Z.-Z. Bai, X. Yang, On semi-convergence of parameterized Uzawa methods for singular saddle point problems, Linear Algebra Appl. 431 (2009) 808–817). In this paper, we discuss the inexact Uzawa method, which covers the Uzawa method, the preconditioned Uzawa method, and the parameterized Uzawa method to solve the singular saddle point problems. We prove the semi-convergence result under restrictions by verifying two necessary and sufficient conditions, that is, all elementary divisors associated with the eigenvalue 1 of its iterative matrix are linear, and the pseudospectral radius of the iterative matrix is less than 1. Sufficient conditions for the semiconvergence of several Uzawa-type methods are also provided. In addition, numerical examples are given to demonstrate the semi-convergence of Uzawa-type methods. © 2013 Elsevier B.V. All rights reserved.
1. Introduction This paper provides the semi-convergence analysis of Uzawa-type methods for singular saddle point problems (S2 P2 ). Consider the iterative solution of a system of linear equations with 2 × 2 block structure,
A B
BT 0
x y
=
f g
,
(1.1)
where A ∈ Rn×n , B ∈ Rm×n , f ∈ Rn and g ∈ Rm , (m ≤ n). Generally speaking, both A and B are large and sparse. We denote the range and the null spaces of A by R(A) and N (A), and the conjugate transpose and transpose of A by A∗ and AT , respectively. Systems of linear equations with the form (1.1) are called saddle point problems (SP2 ), which arise in many scientific and engineering applications [1,2], such as computational fluid dynamics, constrained optimization, incompressible elasticity, circuit analysis, structural analysis, and so forth. Many iterative methods have been proposed in the literature when the linear system (1.1) is invertible. For example, Uzawa-type methods [3–12], matrix splitting iterative methods [13–19], relaxation iterative methods [5,20,21], Krylov subspace iterative methods, (often using block-diagonal, block-tridiagonal, constraint, SOR and HSS preconditioners) [22–24,1,25–32], and iterative null space methods [33,34]. Golub, Wu and Yuan [20] proposed a SOR-like method for
∗
Corresponding author at: School of Mathematical Sciences, Fudan University, Shanghai, 200433, PR China. Tel.: +86 21 65642347. E-mail addresses: [email protected] (N. Zhang), [email protected] (T.-T. Lu), [email protected], [email protected] (Y. Wei).
0377-0427/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cam.2013.05.015
N. Zhang et al. / Journal of Computational and Applied Mathematics 255 (2014) 334–345
335
solving S2 P2 (1.1), and discussed its convergence and choice of optimal parameters. Bai, Parlett and Wang [5] presented a generalized SOR (GSOR) method, which is a two-parameter extension of the SOR-like methods. The Uzawa method is a special case of the GSOR method. It can be found from the asymptotic convergence rate of the Uzawa method that the GSOR method is much more effective than the Uzawa method. Chen and Jiang [35] extended these methods and proposed a class of generalized inexact parameterized iterative schemes for solving saddle point problems. Shao, Shen and Li [21] gave another generalized SOR method, which is also an extension of the SOR-like methods. Bai and Wang [4] further developed the parameterized inexact Uzawa (PIU) method to large sparse generalized saddle point problems. The PIU method is a two-parameter generalization of the inexact Uzawa method, and it is shown [4] that the asymptotic convergence rate of the quasi-optimal PIU method is faster than that of the quasi-optimal inexact Uzawa method. However, both methods have the same computational complexity. In the linear system (1.1), if B is rank-deficient, then the coefficient matrix is singular, and (1.1) is called a singular saddle point problem (S2 P2 ). Recently, several authors [36,26,37–43] have paid attention to the iterative methods for S2 P2 . Zheng, Bai and Yang [43] applied the parameterized Uzawa method to solve S2 P2 . Assume that A is symmetric positive definite, and B is rank-deficient, Zheng, Bai and Yang gave sufficient conditions for the semi-convergence of the parameterized Uzawa method, and determined the optimal iteration parameters. Wu, Silva, and Yuan [44] developed the conjugate gradient method to solve singular saddle point problems, by transforming S2 P2 into a smaller nonsingular SP2 . Cao [37] compared the convergence performance of Krylov subspace methods for solving S2 P2 and the smaller nonsingular SP2 . With a combination of analytic and empirical results, he investigated the effects of singularity and nonsingularity on the convergence, and observed that the convergence behavior of S2 P2 is significantly better than that of the corresponding nonsingular SP2 . Ma and Zhang [38] studied block-diagonally preconditioned parameterized inexact Uzawa methods for S2 P2 , and gave the corresponding semi-convergence analysis. Moreover, Zhang and Shen [42] studied constraint preconditioners for solving singular saddle point problems, and proved that for solving S2 P2 by the preconditioned GMRES methods with constraint preconditioners, GMRES will determine the least squares solutions at breakdown. Zhang and Wei [39–41] investigated general stationary iterative methods and GMRES methods for solving S2 P2 . In this paper we investigate Uzawa-type methods, which cover the Uzawa method, the preconditioned Uzawa method, and the parameterized Uzawa method for solving singular saddle point problems. Since several Uzawa-type methods are special cases of the abstract inexact Uzawa method [6], we mainly concern ourselves with the semi-convergence of the inexact Uzawa method for solving S2 P2 (1.1). We assume the singular saddle point problem (1.1) satisfying the following assumptions in this paper: • Matrix B is rank-deficient, i.e., rank(B) < m; • Matrix A is positive real, i.e., v T Av > 0, for any nonzero v ∈ Rn ; • (1.1) is consistent, i.e., the solution of (1.1) exists. These assumptions are satisfied when the partial differential equations are discretized by finite difference methods or finite element methods [28,44]. The rest of this paper is organized as follows. In Section 2, we consider how to solve the equivalent linear system of (1.1) and introduce Uzawa-type methods. In Section 3, we make a semi-convergence analysis of Uzawa-type methods for solving singular saddle point problems and discuss the optimal parameters. In Section 4, some numerical examples are given to demonstrate the semi-convergence of Uzawa-type methods. 2. Uzawa-type methods To solve (1.1), we can solve the equivalent linear system instead of the original linear system:
BT 0
A −B
BT 0
A −B
x y
=
f
−g
.
(2.1)
Denote
A≡
,
X≡
x , y
F ≡
f
−g
.
Then (2.1) can be rewritten as
AX = F .
(2.2)
It is well known that the coefficient matrix of (1.1) is indefinite. Whether the coefficient matrix of (1.1) is symmetric or not, the coefficient matrix A of (2.1) is nonsymmetric. Although A loses symmetry, similar to the nonsingular case (cf. [27]), it retains some good properties as follows. Lemma 2.1 ([27]). Assume that A is positive real, then (1) A is range-symmetric, that is, R(A) = R(AT ), or equivalently, N (A) = N (AT ). (2) A is semi-positive real: v T Av ≥ 0, for any v ∈ Rn+m . (3) A is positive semi-stable, i.e., the eigenvalues of A have nonnegative real part.
336
N. Zhang et al. / Journal of Computational and Applied Mathematics 255 (2014) 334–345
For the coefficient matrix A of the linear equations (2.1), we utilize the matrix splitting [35]:
A≡
BT 0
A −B
=
A + Q1 −B
0 Q2
−
Q1 0
−BT Q2
≡ M − P,
(2.3)
where A + Q1 is nonsingular and Q2 is symmetric positive definite, then it yields the inexact Uzawa method [6]: xk+1 = xk + (A + Q1 )−1 (f − Axk − BT yk ), yk+1 = yk + Q2−1 (Bxk+1 − g ),
(2.4)
and the iterative matrix T is:
T = M −1 P = Denote Xk ≡
(A + Q1 )−1 Q1 −1 Q2 B(A + Q1 )−1 Q1
−(A + Q1 )−1 BT . I − Q2−1 B(A + Q1 )−1 BT
(2.5)
k x
yk
, then (2.4) can be rewritten as
Xk+1 = T Xk + M−1 F .
(2.6)
With different choices of Q1 and Q2 , (2.4) covers several Uzawa-type methods as follows: • Uzawa method [3]: Q1 = 0, Q2 = τ1 I (τ > 0), xk+1 = xk + A−1 (f − Axk − BT yk ), yk+1 = yk + τ (Bxk+1 − g ).
(2.7)
• Preconditioned Uzawa method [7]: Q1 = 0, k+1 x = xk + A−1 (f − Axk − BT yk ), k+1 y = yk + Q2−1 (Bxk+1 − g ). • Parameterized Uzawa method, or, the generalized successive overrelaxation (GSOR) method [5,43]: Q1 = Q2 = τ1 Q , where 0 < ω ≤ 1, τ > 0, and Q is symmetric positive definite, k+1 x = (1 − ω)xk + ωA−1 (f − BT yk ), yk+1 = yk + τ Q −1 (Bxk+1 − g ).
(2.8) 1
ω
A − A,
(2.9)
3. Semi-convergence analysis In this section, we investigate the convergence analysis of Uzawa-type methods for solving singular saddle point problems. Let us recall some classical results on stationary [40] or nonstationary [45] iterative methods for solving the singular linear equations. For a matrix B ∈ Rn×n , the smallest nonnegative integer j such that rank(Bj ) = rank(Bj+1 ) is called the index of B, and denoted by j = index(B). Thus index(B) is the size of the largest Jordan block corresponding to the zero eigenvalue of B. Moreover, the Drazin inverse [46] of B is defined as the unique matrix BD which satisfies three matrix equations: BD BBD = BD ,
BD B = BBD ,
Bj+1 BD = Bj .
If j = 1, then the Drazin inverse BD reduces to the group inverse, which denoted by Bg . We now introduce the semi-convergence for solving singular linear systems, which has been studied by several authors (cf. [36,47–49]). Definition 3.1 ([5,50,43]). The iterative method (2.6) is semi-convergent if, for any initial vector X0 ∈ Rm+n , the iteration sequence {Xk } produced by the iterative method (2.6) semi-converges to a solution X∗ of singular linear equations (2.2). Moreover, it holds that X∗ = (I − T )g M −1 F + [I − (I − T )g (I − T )]X0 . For the splitting (2.3) and the associated stationary iteration (2.6), it is well known that if A is nonsingular, then the stationary iteration converges if and only if the spectral radius of the iterative matrix T satisfies ρ(T ) < 1, thus limk→∞ T k = 0. As A is singular, then T has eigenvalue 1 and ρ(T ) cannot be smaller than 1. For the iterative matrix T , we introduce its pseudo-spectral radius ν(T ) [47],
ν(T ) = max{|λ| : λ ∈ σ (T ) \ {1}} where σ (T ) is the set of eigenvalues of matrix T . Now we introduce the necessary and sufficient conditions for the semi-convergence of (2.6). Theorem 3.2 ([50]). The iterative method (2.6) is semi-convergent, if and only if index(I − T ) = 1 and ν(T ) < 1.
(3.1)
N. Zhang et al. / Journal of Computational and Applied Mathematics 255 (2014) 334–345
337
Remark 3.3. If ν(T ) < 1 and index(I − T ) = 1 hold, then T is called semi-convergent ([50, p.198]), while limk→∞ T I − (I − T )g (I − T ). Thus the iteration (2.6) is semi-convergent if and only if T is semi-convergent.
k
=
Remark 3.4. It is known that the magnitude of ν(T ) < 1 determines the convergence rate for {T k }, and the smaller ν(T ) is, the higher asymptotic rate of convergence − ln ν(T ) is. For the iteration (2.6), we require the pseudo-spectral radius of its iterative matrix as small as possible. In Section 4, we discuss the optimal parameters which minimize the pseudo-spectral radius ν(T ). 3.1. The conditions for index(I − T ) = 1 In this section we present the conditions for index(I − T ) = 1. Lemma 3.5 ([40]). Index(I − T ) = 1 holds if and only if, for any 0 ̸= Y ∈ R(A), Y ̸∈ N (AM −1 ). Proof. (⇒) For any 0 ̸= Y ∈ R(A), AM −1 Y ̸= 0, we prove index(I − T ) = 1 by contradiction. Suppose index(I − T ) > 1, then there exists a nonzero X ∈ Rm+n such that M −1 AX ̸= 0 and (M −1 A)2 X = 0. Let Y = AX , then Y ̸= 0, (M−1 A)2 X = M−1 AM−1 Y = 0 ⇒ AM−1 Y = 0, which is in contradiction with 0 ̸= Y ∈ R(A), AM−1 Y ̸= 0. (⇐) If index(I − T ) = 1. We can show that for any 0 ̸= Y ∈ R(A), AM−1 Y ̸= 0. Suppose that there exists a Y = AX ̸= 0 such that AM −1 Y = 0. AX ̸= 0 ⇒ M −1 AX ̸= 0. (M −1 A)2 X = M −1 (AM −1 Y ) = 0, then rank[(M −1 A)2 ] < rank(M −1 A), which contradicts with index(I − T ) = 1. Theorem 3.6. Assume that A is positive real, A + Q1 is nonsingular, and Q2 is symmetric positive definite, then index(I − T ) = 1. Proof. Let X = (ξ T , ηT )T , where ξ ∈ Rn , η ∈ Rm such that
Y = AX =
Aξ + BT η −Bξ
̸= 0.
Notice that
M=
A + Q1 −B
0 Q2
,
M
−1
(A + Q1 )−1 = −1 Q2 B(A + Q1 )−1
0 Q2−1
,
and
AM−1 =
A(A + Q1 )−1 + BT Q2−1 B(A + Q1 )−1 −B(A + Q1 )−1
BT Q2−1 0
,
then
(A + BT Q2−1 B)(A + Q1 )−1 (Aξ + BT η) − BT Q2−1 Bξ AM Y = . −B(A + Q1 )−1 (Aξ + BT η) −1
Denote φ = (A + Q1 )−1 (Aξ + BT η), then
AM−1 Y =
(A + BT Q2−1 B)φ − BT Q2−1 Bξ . −Bφ
To prove index(I − T ) = 1, it is sufficient to prove Y ̸∈ N (AM −1 ) by Lemma 3.5. We fulfill the proof based on the following three cases. 0
Case (a) Aξ + BT η = 0: Y = −Bξ
̸= 0, then Bξ ̸= 0. AM−1 Y =
so is Q2 , then it is easy to see that Y ̸∈ N (AM −1
−1
−BT Q2−1 Bξ 0
, note that Q2 is symmetric positive definite,
).
Case (b) Aξ + B η ̸= 0, Bφ ̸= 0: In this case, it is obvious that Y ̸∈ N (AM −1 ) holds. Case (c) Aξ + BT η ̸= 0, Bφ = 0: T
AM−1 Y =
Aφ − BT Q2−1 Bξ 0
.
To prove Y ̸∈ N (AM −1 ), it is sufficient to prove Aφ − BT Q2−1 Bξ ̸= 0. Next we provide the proof for Aφ − BT Q2−1 Bξ ̸= 0 by contradiction. Assume Aφ − BT Q2−1 Bξ = 0, i.e., φ = A−1 BT Q2−1 Bξ . It follows from Aξ + BT η ̸= 0 that φ ̸= 0, which leads to Bξ ̸∈ N (A−1 BT Q2−1 ).
(3.2)
On the other hand, since Bφ = BA−1 BT Q2 Bξ = 0, we have −1
Bξ ∈ N (BA−1 BT Q2−1 ).
(3.3)
338
N. Zhang et al. / Journal of Computational and Applied Mathematics 255 (2014) 334–345
It follows from (3.2) and (3.3) that rank(A−1 BT Q2−1 ) > rank(BA−1 BT Q2−1 ).
(3.4)
In fact, rank(A−1 BT Q2−1 ) = rank(BT ) = rank(B), and notice that A is positive real, then rank(BA−1 BT Q2−1 ) = rank(BA−1 BT ) = rank(B). The above two equations yield rank(A−1 BT Q2−1 ) = rank(BA−1 BT Q2−1 ), which is in contradiction with (3.4).
3.2. The conditions for ν(T ) < 1
u
Let T
v
=λ
u u n m v , v ̸= 0, where u ∈ C , v ∈ C , that is,
Q1 u − BT v = λ(A + Q1 )u, Q2 v = −λBu + λQ2 v.
(3.5)
λ From (3.5), when λ ̸= 1, it holds that v = λ− Q −1 Bu and 1 2
(Q1 − λA − λQ1 )u =
λ T −1 B Q2 Bu. λ−1
(3.6)
Denote
α=
u∗ Q1 u u∗ u
,
β=
u∗ Au u∗ u
,
γ =
u∗ BT Q2−1 Bu
Multiplying both sides of (3.6) from the left with
λ2 +
u∗ u u∗ u∗ u
.
(3.7)
yields
α γ − 2α − β λ+ = 0. α+β α+β
(3.8)
Lemma 3.7. Suppose that A is positive real, A + Q1 is nonsingular, and Q2 is symmetric positive definite, then u = 0, if and only if λ = 1. Proof. If u = 0, then it follows from (3.5) that (1 −λ)Q2 v = 0, notice that v ̸= 0, or else (uT , v T )T = 0, so λ = 1. Conversely, if λ = 1, then −BT v = Au and Bu = 0. Thus B(−A−1 BT v) = Bu = 0, i.e., (BT v)∗ A−1 (BT v) = 0. Notice that A is positive real, so is A−1 , then BT v = 0, and Au = 0, thus u = 0. Lemma 3.8. Assume that A is positive real, Q1 is symmetric positive semi-definite, and Q2 is symmetric positive definite. If u ∈ N (B), then |λ| < 1. Proof. From (3.5), it holds that (1 − λ)Q2 v = 0, and it follows from Lemma 3.7 that λ ̸= 1, then v = 0. Again from (3.5) α we have (1 − λ)Q1 u = λAu, which yields (α + β)λ = α . Notice that α ≥ 0, and the real part of β is positive, so |λ| = |α+β| < 1. To discuss the conditions for |λ| < 1 with u ̸∈ N (B), we introduce a result about the evaluation of the roots of a quadratic equation as follows. Lemma 3.9 ([51]). Both roots of the quadratic equation x2 −px+q = 0 are less than 1 in modulus if and only if |p−¯pq| < 1−|q|2 , where p¯ is the conjugate number of p. If p and q are real, then the necessary and sufficient conditions for both roots of the quadratic equation x2 − px + q = 0 are less than 1 in modulus is |q| < 1 and |p| < 1 + q [19]. Lemma 3.10. Assume A is positive real, Q1 is symmetric positive semi-definite, and Q2 is symmetric positive definite. If u ̸∈ N (B), then |λ| < 1, if and only if
γ < 2c +
4α c 2 c 2 + d2
,
where β = c + id, or, c and d are the real part and imaginary part of β , respectively, and i =
(3.9)
√ −1.
N. Zhang et al. / Journal of Computational and Applied Mathematics 255 (2014) 334–345
339
Proof. Since u ̸= 0, from Lemma 3.7 it holds that λ ̸= 1, which yields Eq. (3.8)), i.e.,
λ2 +
γ − 2α − β α λ+ = 0. α+β α+β
From Lemma 3.9, we have
γ − 2α − β α 2 γ − 2α − β¯ α . − |λ| < 1 ⇔ <1− α+β α +β α +β α + β¯ By a little algebra, we obtain
|λ| < 1 ⇔ (γ − α)β¯ − βα − β β¯ < |α + β|2 − α 2 , i.e.,
|λ| < 1 ⇔ (α − γ )(c − id) + α(c + id) + c 2 + d2 < |α + c + id|2 − α 2 . Denote s = c 2 + d2 + 2α c, and by some simple algebra operations, we obtain
|λ| < 1 ⇔ |s − γ c + iγ d| < s ⇔ |s − γ c + iγ d|2 < s2 ⇔γ 2 c 2 + γ 2 d2 − 2γ cs < 0. Note that Bu ̸= 0, so γ > 0, then
|λ| < 1 ⇔ γ c 2 + γ d2 − 2cs < 0 ⇔γ (c 2 + d2 ) < 2c (c 2 + d2 + 2α c ) 4α c 2 , ⇔γ < 2c + 2 c + d2 which finishes the proof.
Together with Lemmas 3.7, 3.8 and 3.10, we prove the following result. Theorem 3.11. Assume that A is positive real, Q1 is symmetric positive semi-definite, and Q2 is symmetric positive definite. Then ν(T ) < 1 if and only if, when u ̸∈ N (B), (3.9) holds. Corollary 3.12. Suppose that A and Q2 is symmetric positive definite, and Q1 is symmetric positive semi-definite. Then ν(T ) < 1 if and only if, when u ̸∈ N (B),
γ < 2c + 4α.
(3.10)
3.3. Sufficient conditions for the semi-convergence of Uzawa-type methods The following result gives sufficient conditions for the semi-convergence of the inexact Uzawa method. Theorem 3.13. Assume that A is positive real, Q1 is symmetric positive semi-definite, and Q2 is symmetric positive definite. If
λmax (BT Q2−1 B) < λmin (A + AT ) +
λmin (Q1 )λ2min (A + AT ) , ∥A∥22
(3.11)
then the inexact Uzawa method (2.6) is semi-convergent, where λmax and λmin denote the largest and smallest eigenvalues of the corresponding matrix, respectively. Proof. By Theorems 3.6 and 3.11, it only needs to prove (3.9) when u ̸∈ N (B). Notice that
γ =
u∗ (BT Q2−1 B)u u∗ u
≤ λmax (BT Q2−1 B),
then it follows from (3.11) that
γ < λmin (A + AT ) +
λmin (Q1 )λ2min (A + AT ) . ∥A∥22
It is obvious that
β = c + id =
u∗ Au u∗ u
,
i.e., c 2 + d2 = |β|2 ≤ ∥A∥22 ,
(3.12)
340
N. Zhang et al. / Journal of Computational and Applied Mathematics 255 (2014) 334–345
and c=
1 u∗ (A + AT )u u∗ u
2
≥
1 2
λmin (A + AT ).
Together with the above two inequalities and α =
u∗ Q1 u u∗ u
≥ λmin (Q1 ), we have
4α c 2 λmin (Q1 )λ2min (A + AT ) ≤ 2c + . c 2 + d2 ∥A∥22
λmin (A + AT ) +
4α c 2 , c 2 +d2
From (3.12) and (3.13), it is clear that γ < 2c +
which completes the proof.
(3.13)
Corollary 3.14. Suppose that A and Q2 are symmetric positive definite, Q1 is symmetric positive semi-definite, and
λmax (BT Q2−1 B) < 2λmin (A) + 4λmin (Q1 )
λ2min (A) . λ2max (A)
Then the inexact Uzawa method (2.6) is semi-convergent. The following results give sufficient conditions for the semi-convergence of the Uzawa method. Theorem 3.15. Assume that A is positive real, then the Uzawa method (2.7) is semi-convergent, if 0<τ <
1
. λmax (A + AT )−1 BT B
(3.14)
Proof. By Theorems 3.6 and 3.11, we only need to prove (3.9) when u ̸∈ N (B). Notice that
α = 0,
β=
u∗ Au u∗ u
,
γ =τ
u∗ BT Bu u∗ u
,
then (3.9) ⇔ τ
u∗ BT Bu u∗ u
<
u∗ (A + AT )u
u∗ u ⇔ u (τ B B − (A+ AT ))u < 0 ∗
T
1 1 1 1 ⇔ u∗ (A + AT ) 2 I − τ (A + AT )− 2 BT B(A + AT )− 2 (A + AT ) 2 u > 0. 1
1
The above inequality holds if τ λ[(A + AT )− 2 BT B(A + AT )− 2 ] < 1. In fact, 1
1
λ[(A + AT )− 2 BT B(A + AT )− 2 ] = λ[(A + AT )−1 BT B] ≤ λmax [(A + AT )−1 BT B], then from (3.14) it is easy to see that (3.9) holds.
Corollary 3.16. Assume that A is symmetric positive definite, then the Uzawa method (2.7) is semi-convergent, if 0<τ <
2
λmax (A−1 BT B)
.
(3.15)
The following results give sufficient conditions for the semi-convergence of the preconditioned Uzawa method. Theorem 3.17. Assume that A is positive real, Q2 is symmetric positive definite, then the preconditioned Uzawa method (2.8) is semi-convergent, if
λmax [(A + AT )−1 BT Q2−1 B] < 1. Proof. The proof is similar to that of Theorem 3.15.
(3.16)
Corollary 3.18. Suppose that A and Q2 are symmetric positive definite, then the preconditioned Uzawa method (2.8) is semiconvergent, if
λmax (A−1 BT Q2−1 B) < 2.
(3.17)
Theorem 3.19 gives sufficient conditions for the semi-convergence of the parameterized Uzawa method, which is presented in [43], we give another proof for completeness.
N. Zhang et al. / Journal of Computational and Applied Mathematics 255 (2014) 334–345
341
Theorem 3.19 ([43]). Assume that A and Q are symmetric positive definite, then the parameterized Uzawa method (2.9) is semiconvergent, if 0 < ω < 2,
0<τ <
2(2 − ω)
.
ωλmax (Q −1 BA−1 BT )
(3.18)
Lemma 3.20. Suppose that A and Q are symmetric positive definite, then for the parameterized Uzawa method (2.9), it holds that ν(T ) < 1, if 0 < ω < 2,
γ < 2c + 4α.
(3.19)
Proof. Recall Eq. (3.7), and c is the real part of β ,
α=
1 − ω u∗ Au
ω
u∗ u
,
β=c=
u∗ Au u∗ u
,
γ =τ
u∗ BT Q −1 Bu u∗ u
.
If u = 0, then from Lemma 3.7 we obtain λ = 1.
α If u ̸= 0, and u ∈ N (B), then from the proof of Lemma 3.8 we obtain λ = α+ = 1 − ω. If 0 < ω < 2, then it holds that c |λ| < 1.
If u ̸∈ N (B), taking the same proof of Lemma 3.10, then we obtain |λ| < 1 from γ < 2c + 4α , which finishes the proof.
Proof of Theorem 3.19. By Theorem 3.6 and Lemma 3.20, we only need to prove γ < 2c + 4α when u ̸∈ N (B). Denote 2(2−ω) s = ω , then
γ < 2c + 4α ⇔ u∗ (sA − τ BT Q −1 B)u > 0 1 1 1 1 ⇔ u∗ A 2 (sI − τ A− 2 BT Q −1 BA− 2 )A 2 u > 0. 1
1
1
1
The above inequality holds if τ λmax (A− 2 BT Q −1 BA− 2 ) < s, and λmax (A− 2 BT Q −1 BA− 2 ) = λmax (Q −1 BA−1 BT ), then from (3.18), it is easy to see that γ < 2c + 4α holds, which completes the proof. 3.4. The optimal parameters of Uzawa-type methods Among the Uzawa-type methods we considered, both the Uzawa method (2.7) and the parameterized Uzawa method (2.9) have parameters to be chosen. We now discuss the optimal parameters which minimize the pseudo-spectral radius ν(T ). For the parameterized Uzawa method, Zheng, Bai and Yang [43] gave the optimal parameters ωopt and τopt , as shown by the following theorem. Theorem 3.21 ([43]). Let A and Q be symmetric positive definite. Denote the smallest and the largest nonzero eigenvalues of the matrix Q −1 BA−1 BT by λmin and λmax , respectively. Then for the parameterized Uzawa method (2.9), the optimal parameters ωopt and τopt which minimize ν(T ) are
√ 4 λmin λmax , ωopt = √ √ ( λmin + λmax )2
1
τopt = √ , λmin λmax
(3.20)
and
√ √ λmax − λmin min ν(T ) = √ . √ λmax + λmin
(3.21)
The following theorem presents the optimal parameter τopt for the Uzawa method (2.7) to solve S2 P2 (1.1), and the result is similar to that of the nonsingular SP2 . Theorem 3.22. Assume A is symmetric positive definite. Denote the smallest and the largest nonzero eigenvalues of the matrix A−1 BT B by λmin and λmax , respectively. Then for the Uzawa method (2.7), the optimal parameter τopt which minimizes ν(T ) is
τopt =
2
,
(3.22)
λmax − λmin . λmax + λmin
(3.23)
λmin + λmax
and min ν(T ) =
342
N. Zhang et al. / Journal of Computational and Applied Mathematics 255 (2014) 334–345
Proof. Lemma 3.7 implies u = 0 if and only if λ = 1, which has nothing to do with ν(T ). Assume u ̸= 0 and recall (3.7) and (3.8), we have
α = 0,
β=c=
u∗ Au u∗ u
,
γ =τ
u∗ BT Bu u∗ u
,
and
λ2 +
γ −c c
λ = 0.
Then λ = 0 or λ = 1 − Denote θ =
u∗ BT Bu , u∗ Au
γ
. Notice that γ ̸= 0, or else it yields λ = 1, so Bu ̸= 0.
c
then λ = 1 − τ θ . Since 1
1
1
1
uˆ ∗ (A− 2 BT BA− 2 )ˆu (A− 2 uˆ )∗ BT B(A− 2 uˆ ) = , 0<θ = = 1 1 uˆ ∗ uˆ uˆ ∗ uˆ (A 2 u)∗ (A 2 u) u∗ BT Bu
1
1
1
where uˆ = A 2 u, and notice λ(A− 2 BT BA− 2 ) = λ(A−1 BT B), then λmin ≤ θ ≤ λmax . Now, for the following ‘‘min–max’’ problem min τ
max
λmin ≤θ≤λmax
|λ|,
it is easy to see that the optimal parameter τopt is attained when 1 − τ λmax = −(1 − τ λmin ), that is,
τopt =
2
λmin + λmax
,
and
min ν(T ) = 1 − τopt λmin =
λmax − λmin . λmax + λmin
Remark 3.23. In this subsection we discuss the optimal parameters which minimize the pseudo-spectral radius ν(T ). Notice that if Q is the identity matrix, then λmin and λmax are the same as Theorems 3.21 and 3.22, and it is easy to check that
√ √ λmax − λmin λmax − λmin . < √ √ λ λmax + λmin max + λmin
This means min ν(T ) of Theorem 3.21 is smaller than that of Theorem 3.22, or the parameterized Uzawa method is faster than the Uzawa method. 4. Numerical experiments In this section we compare the performance of several Uzawa-type methods. All the computations are implemented in MATLAB on a PC computer with Genuine Intel (R) CPU 1.83 GHz, and 512 MB memory. Consider the Navier–Stokes equations,
−ν △ u + (u · ▽)u + ▽p = f −▽·u=0
in Ω , in Ω ,
(4.1)
where Ω is an open bounded domain in R2 or R3 , the vector field u represents the velocity in Ω , p denotes pressure, and the scalar ν is the viscosity, which is inversely proportional to the Reynolds number. Linearizing Eq. (4.1) by Picard iteration, we get the Oseen equations,
−ν △ u + (w · ▽)u + ▽p = f −▽·u=0
in Ω , in Ω ,
(4.2)
where w satisfies such that ▽ · w = 0. The test problem is a ‘‘leaky’’ two-dimensional lid-driven cavity problem in the square domain: Ω = (0 < x < 1 : 0 < y < 1). Let u = (u, v)T denote the velocity field and w = (a, b)T the wind. The boundary conditions are u = v = 0 on the three fixed walls (x = 0, y = 0, x = 1), and u = 1, v = 0 on the moving wall (y = 1). We take constant ‘‘wind’’ a = 1, b = 2. To discretize (4.2), we use the ‘‘marker and cell’’ (MAC) finite difference scheme [28,52]. Divide Ω into a uniform n × n grid of cells of width h = 1/n. The discrete velocities and pressures are defined on a staggered grid in which the discrete values of u lie at the centers of the cell boundaries orthogonal to the x-axis, the discrete values of v lie at the centers of the cell boundaries orthogonal to the y-axis, and the discrete pressures lie at the cell centers. We obtain the matrix representation of the Oseen equations (4.2),
AX ≡
A B
BT 0
u p
=
f , 0
(4.3)
N. Zhang et al. / Journal of Computational and Applied Mathematics 255 (2014) 334–345
343
Table 1 CPU time, iteration numbers, and residual with θ = 0.01. Iterative method
CPU (s)
IT
RES
Preconditioned GMRES Uzawa Preconditioned Uzawa Parameterized Uzawa
11.43
109 Divergent 399 Divergent
9.13e-09
15.91
9.88e-09
Table 2 CPU time, iteration numbers, and residual with θ = 0.001. Iterative method
CPU (s)
IT
RES
Preconditioned GMRES Uzawa Preconditioned Uzawa Parameterized Uzawa
9.30 8.63 13.64 61.65
38 72 275 500
6.90e-09 9.42e-09 9.90e-09 3.44e-03
where
A=
F1 0
0 F2
∈ R2n(n−1)×2n(n−1) ,
2 B = (B1 , B2 ) ∈ Rn ×2n(n−1) ,
Fi = ν Ai + Ni ∈ Rn(n−1)×n(n−1) ,
(i = 1, 2), u ∈ R2n(n−1) , p ∈ Rn . For convenience we denote l = 2n(n − 1), m = n2 in this section. The coefficient matrix of linear systems (4.3) satisfies the 2
following properties: A is nonsymmetric and positive real, rank(B) = m−1, rank(A) = l+m−1, N (BT ) = span{em }, em = (1, 1, . . . , 1)T ∈ Rm . To ensure the (1, 1)-block matrix is symmetric positive definite, finally we take the test coefficient matrix as follows,
A¯ =
BT 0
A¯ −B
=
1
(A + AT ) 2 −B
BT 0
.
(4.4)
We compare Uzawa-type methods with the preconditioned GMRES method, and the preconditioner is the block triangular preconditioner [28],
¯
A
Q=
0
BT 1 2 h I
.
(4.5)
ν
For the preconditioned Uzawa method (2.8) and the parameterized Uzawa method (2.9), Q2 and Q are taken as the identity matrix I. Let n = 32, ν = 0.5, then the smallest and the largest nonzero eigenvalues of the matrix A¯ −1 BT B are about 0.0151 and 0.0625, respectively. We take the optimal parameters by Theorems 3.21 and 3.22. For all iterative methods, the initial vector is the zero vector. When the iteration number is k, denote the k-th iterative residual by r k , the stopping criteria is:
∥r k ∥2 ≤ 10−8 . ∥r 0 ∥2 We use the preconditioned restarted GMRES(20) in our numerical tests. The iterative number k of GMRES means k = (i − 1) × 20 + j, where i is the number of restarting, and j is the iteration number of the last restarting, respectively. Notice that it is necessary to compute A¯ −1 for all Uzawa-type methods, and for the preconditioned GMRES, to compute the matrix–vector products with Q−1 , it is still necessary to compute A¯ −1 . Thus, for all iterative methods, we use the MATLAB code cholinc (A¯ , θ ) to approximate A¯ −1 , where θ is the dropping tolerance. With respect to different dropping tolerance θ , we list the CPU time (denoted as CPU), iteration numbers (denoted as IT), and the norm of the residual (denoted as RES) for the above different iterative methods in the tables. Remark 4.1. It can be seen from Tables 1–4, for the Uzawa method, especially for the parameterized Uzawa method, the computation results depend on the exact computation of A¯ −1 , and for the preconditioned GMRES and the preconditioned Uzawa, the computational results do not rely on the exact computation of A¯ −1 . Moreover, the preconditioned GMRES and the preconditioned Uzawa do not require the estimation of the smallest and the largest nonzero eigenvalues of the matrix A¯ −1 BT B. Acknowledgments We would like to express our sincere thanks to our editor Prof. Michael K. Ng and three anonymous reviewers for their valuable suggestions that greatly improved the representation of this paper. Finally, we thank Prof. Eric King-wah Chu for carefully reading our draft and his many valuable suggestions.
344
N. Zhang et al. / Journal of Computational and Applied Mathematics 255 (2014) 334–345 Table 3 CPU time, iteration numbers, and residual with θ = 0.0001. Iterative method
CPU (s)
IT
RES
Preconditioned GMRES Uzawa Preconditioned Uzawa Parameterized Uzawa
8.75 7.67 10.90 61.41
26 22 153 500
1.92e-09 8.96e-09 9.96e-09 6.99e-04
Table 4 CPU time, iteration numbers, and residual with θ = 0.00001. Iterative method
CPU (s)
IT
RES
Preconditioned GMRES Uzawa Preconditioned Uzawa Parameterized Uzawa
8.72 7.64 7.84 61.55
24 16 23 500
5.54e-09 8.89e-09 9.75e-09 1.49e-05
The first author was supported by Zhejiang Provincial Natural Science Foundation of China under grant No. Y1110451 and National Natural Science Foundation of China under grant No. 61002039. The third author was supported by the National Natural Science Foundation of China under grant 11271084, the Doctoral Program of the Ministry of Education under grant 20090071110003, 973 Program Project (No. 2010CB327900) and Shanghai Education Committee. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38]
Z.-Z. Bai, Structured preconditioners for nonsingular matrices of block two-by-two structures, Math. Comp. 75 (2006) 791–815. M. Benzi, G.H. Golub, J. Liesen, Numerical solution of saddle point problems, Acta Numer. 14 (2005) 1–137. K. Arrow, L. Hurwicz, H. Uzawa, Studies in Nonlinear Programming, Stanford University Press, Stanford, CA, 1958. Z.-Z. Bai, Z.-Q. Wang, On parameterized inexact Uzawa methods for generalized saddle point problems, Linear Algebra Appl. 428 (2008) 2900–2932. Z.-Z. Bai, B.N. Parlett, Z.-Q. Wang, On generalized successive overrelaxation methods for augmented linear systems, Numer. Math. 102 (2005) 1–38. 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. H.C. Elman, G.H. Golub, Inexact and preconditioned Uzawa algorithms for saddle point problems, SIAM J. Numer. Anal. 31 (1994) 1645–1661. Q. Hu, J. Zou, An iterative method with variable relaxation parameters for saddle point problems, SIAM J. Matrix Anal. Appl. 23 (2001) 317–338. Q. Hu, J. Zou, Two new variants of nonlinear inexact Uzawa algorithms for saddle point problems, Numer. Math. 93 (2002) 333–359. Y. Lin, Y. Wei, Fast corrected Uzawa methods for solving symmetric saddle point problems, Calcolo 43 (2006) 65–82. J. Lu, Z. Zhang, A modified nonlinear inexact Uzawa algorithm with a variable relaxation parameter for the stabilized saddle point problem, SIAM J. Matrix Anal. Appl. 31 (2010) 1934–1957. W. Zulehner, Analysis of iterative methods for saddle point problems: a unified approach, Math. Comp. 71 (2002) 479–505. 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, 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. 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. Z.-Z. Bai, Optimal parameters in the HSS-like methods for saddle point problems, Numer. Linear Algebra Appl. 16 (2009) 447–479. M.-Q. Jiang, Y. Cao, On local Hermitian and skew-Hermitian splitting iteration methods for generalized saddle point problems, J. Comput. Appl. Math. 231 (2009) 973–982. X. Peng, W. Li, The alternating-direction iterative method for saddle point problems, Appl. Math. Comput. 216 (2010) 1845–1858. D.M. Young, Iterative Solution for Large Linear Systems, Academic Press, New York, 1971. G.H. Golub, X. Wu, J.-Y. Yuan, SOR-like methods for augmented systems, BIT 41 (2001) 71–85. X. Shao, H. Shen, C. Li, The generalized SOR-like method for the augmented systems, Int. J. Inf. Syst. Sci. 2 (1) (2006) 92–98. 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, 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. Z.-Z. Bai, M.K. Ng, On inexact preconditioners for nonsymmetric matrices, SIAM J. Sci. Comput. 26 (2005) 1710–1724. Z.-Z. Bai, M.K. Ng, Z.-Q. Wang, Constraint preconditioners for symmetric indefinite matrices, SIAM J. Matrix Anal. Appl. 31 (2009) 410–433. Z.-Z. Bai, G.H. Golub, C.-K. Li, Convergence properties of preconditioned Hermitian and skew-Hermitian splitting methods for non-Hermitian positive semidefinite matrices, Math. Comp. 76 (2007) 287–298. M. Benzi, G.H. Golub, A preconditioner for generalized saddle point problems, SIAM J. Matrix Anal. Appl. 26 (2004) 20–41. H.C. Elman, Preconditioning for the steady-state Navier–Stokes equations with low viscosity, SIAM J. Sci. Comput. 20 (4) (1999) 1299–1316. C. Keller, N.I.M. Gould, A.J. Wathen, Constraint preconditioning for indefinite linear systems, SIAM J. Matrix Anal. Appl. 21 (2000) 1300–1317. Y. Lin, Y. Wei, A note on constraint preconditioners for nonsymmetric saddle point problems, Numer. Linear Algebra Appl. 14 (2007) 659–664. J.-Y. Pan, M.K. Ng, Z.-Z. Bai, New preconditioners for saddle point problems, Appl. Math. Comput. 172 (2006) 762–771. E.D. Sturier, J. Liesen, Block-diagonal preconditioners for indefinite linear algebraic systems. Part I: theory, SIAM J. Sci. Comput. 26 (5) (2005) 1598–1619. N.I.M. Gould, M.E. Hribar, J. Nocedal, On the solution of equality constrained quadratic programming problems arising in optimization, SIAM J. Sci. Comput. 23 (2001) 1376–1395. V. Sarin, A. Sameh, An efficient iterative method for the generalized Stokes problem, SIAM J. Sci. Comput. 19 (1998) 206–226. F. Chen, Y.-L. Jiang, A generalization of the inexact parameterized Uzawa methods for saddle point problems, Appl. Math. Comput. 206 (2008) 765–771. Z.-Z. Bai, On semi-convergence of Hermitian and skew-Hermitian splitting methods for singular linear systems, Computing 89 (2010) 171–197. Z. Cao, Comparison of performance of iterative methods for singular and nonsingular saddle point linear systems arising from Navier–Stokes equations, Appl. Math. Comput. 174 (2006) 630–642. H. Ma, N. Zhang, A note on block-diagonally preconditioned PIU methods for singular saddle point problems, Int. J. Comput. Math. 88 (2011) 3448–3457.
N. Zhang et al. / Journal of Computational and Applied Mathematics 255 (2014) 334–345
345
[39] N. Zhang, Y. Wei, Solving EP singular linear systems, Int. J. Comput. Math. 81 (2004) 1395–1405. [40] N. Zhang, Y. Wei, On the convergence of general stationary iterative methods for Range–Hermitian singular linear systems, Numer. Linear Algebra Appl. 17 (2010) 139–154. [41] N. Zhang, A note on preconditioned GMRES for solving singular linear equations, BIT 50 (2010) 207–220. [42] N. Zhang, P. Shen, Constraint preconditioners for solving singular saddle point problems, J. Comput. Appl. Math. 238 (2013) 116–125. [43] B. Zheng, Z.-Z. Bai, X. Yang, On semi-convergence of parameterized Uzawa methods for singular saddle point problems, Linear Algebra Appl. 431 (2009) 808–817. [44] X. Wu, B.P.B. Silva, J.Y. Yuan, Conjugate gradient method for rank deficient saddle point problems, Numer. Algorithms 35 (2004) 139–154. [45] X. Shi, Y. Wei, W. Zhang, Convergence of general nonstationary iterative methods for solving singular linear equations, SIAM J. Matrix Anal. Appl. 32 (2011) 72–89. [46] S. Campbell, C. Meyer, Generalized Inverses of Linear Transformations, Pitman, London, 1979, Dover, New York, 1991; SIAM, Philadelphia, 2008. [47] Z. Cao, On the convergence of iterative methods for solving singular linear systems, J. Comput. Appl. Math. 145 (2002) 1–9. [48] W. Li, Y. Liu, X. Peng, The generalized HSS method for solving singular linear systems, J. Comput. Appl. Math. 236 (2012) 2338–2353. [49] Y. Song, L. Wang, Semiconvergence of P-regular splittings for solving singular linear systems, Calcolo 45 (2008) 247–261. [50] A. Berman, R. Plemmons, Nonnegative Matrices in Mathematical Science, Academic Press, New York, 1979, SIAM, 1994. [51] 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. [52] F. Harlow, J. Welch, Numerical calculation of times-dependent viscous incompressible flow of fluid with free surface, Phys. Fluids 8 (1965) 2182–2189.
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
Semi-convergence analysis of Uzawa methods for singular saddle point problems Naimin Zhang a , Tzon-Tzer Lu b , Yimin Wei c,d,∗ a
School of Mathematics and Information Science, Wenzhou University, Wenzhou, 325035, PR China
b
Department of Applied Mathematics, National Sun Yat-sen University, Kaohsiung, 80424, Taiwan
c
School of Mathematical Sciences, Fudan University, Shanghai, 200433, PR China
d
Shanghai Key Laboratory of Contemporary Applied Mathematics, PR China
article
abstract
info
Article history: Received 19 January 2012 Received in revised form 28 November 2012 MSC: 15A09 65F10 Keywords: Singular linear systems Saddle point problems Uzawa method Semi-convergence
Recently, Zheng, Bai and Yang studied the parameterized Uzawa method for solving singular saddle point problems (B. Zheng, Z.-Z. Bai, X. Yang, On semi-convergence of parameterized Uzawa methods for singular saddle point problems, Linear Algebra Appl. 431 (2009) 808–817). In this paper, we discuss the inexact Uzawa method, which covers the Uzawa method, the preconditioned Uzawa method, and the parameterized Uzawa method to solve the singular saddle point problems. We prove the semi-convergence result under restrictions by verifying two necessary and sufficient conditions, that is, all elementary divisors associated with the eigenvalue 1 of its iterative matrix are linear, and the pseudospectral radius of the iterative matrix is less than 1. Sufficient conditions for the semiconvergence of several Uzawa-type methods are also provided. In addition, numerical examples are given to demonstrate the semi-convergence of Uzawa-type methods. © 2013 Elsevier B.V. All rights reserved.
1. Introduction This paper provides the semi-convergence analysis of Uzawa-type methods for singular saddle point problems (S2 P2 ). Consider the iterative solution of a system of linear equations with 2 × 2 block structure,
A B
BT 0
x y
=
f g
,
(1.1)
where A ∈ Rn×n , B ∈ Rm×n , f ∈ Rn and g ∈ Rm , (m ≤ n). Generally speaking, both A and B are large and sparse. We denote the range and the null spaces of A by R(A) and N (A), and the conjugate transpose and transpose of A by A∗ and AT , respectively. Systems of linear equations with the form (1.1) are called saddle point problems (SP2 ), which arise in many scientific and engineering applications [1,2], such as computational fluid dynamics, constrained optimization, incompressible elasticity, circuit analysis, structural analysis, and so forth. Many iterative methods have been proposed in the literature when the linear system (1.1) is invertible. For example, Uzawa-type methods [3–12], matrix splitting iterative methods [13–19], relaxation iterative methods [5,20,21], Krylov subspace iterative methods, (often using block-diagonal, block-tridiagonal, constraint, SOR and HSS preconditioners) [22–24,1,25–32], and iterative null space methods [33,34]. Golub, Wu and Yuan [20] proposed a SOR-like method for
∗
Corresponding author at: School of Mathematical Sciences, Fudan University, Shanghai, 200433, PR China. Tel.: +86 21 65642347. E-mail addresses: [email protected] (N. Zhang), [email protected] (T.-T. Lu), [email protected], [email protected] (Y. Wei).
0377-0427/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cam.2013.05.015
N. Zhang et al. / Journal of Computational and Applied Mathematics 255 (2014) 334–345
335
solving S2 P2 (1.1), and discussed its convergence and choice of optimal parameters. Bai, Parlett and Wang [5] presented a generalized SOR (GSOR) method, which is a two-parameter extension of the SOR-like methods. The Uzawa method is a special case of the GSOR method. It can be found from the asymptotic convergence rate of the Uzawa method that the GSOR method is much more effective than the Uzawa method. Chen and Jiang [35] extended these methods and proposed a class of generalized inexact parameterized iterative schemes for solving saddle point problems. Shao, Shen and Li [21] gave another generalized SOR method, which is also an extension of the SOR-like methods. Bai and Wang [4] further developed the parameterized inexact Uzawa (PIU) method to large sparse generalized saddle point problems. The PIU method is a two-parameter generalization of the inexact Uzawa method, and it is shown [4] that the asymptotic convergence rate of the quasi-optimal PIU method is faster than that of the quasi-optimal inexact Uzawa method. However, both methods have the same computational complexity. In the linear system (1.1), if B is rank-deficient, then the coefficient matrix is singular, and (1.1) is called a singular saddle point problem (S2 P2 ). Recently, several authors [36,26,37–43] have paid attention to the iterative methods for S2 P2 . Zheng, Bai and Yang [43] applied the parameterized Uzawa method to solve S2 P2 . Assume that A is symmetric positive definite, and B is rank-deficient, Zheng, Bai and Yang gave sufficient conditions for the semi-convergence of the parameterized Uzawa method, and determined the optimal iteration parameters. Wu, Silva, and Yuan [44] developed the conjugate gradient method to solve singular saddle point problems, by transforming S2 P2 into a smaller nonsingular SP2 . Cao [37] compared the convergence performance of Krylov subspace methods for solving S2 P2 and the smaller nonsingular SP2 . With a combination of analytic and empirical results, he investigated the effects of singularity and nonsingularity on the convergence, and observed that the convergence behavior of S2 P2 is significantly better than that of the corresponding nonsingular SP2 . Ma and Zhang [38] studied block-diagonally preconditioned parameterized inexact Uzawa methods for S2 P2 , and gave the corresponding semi-convergence analysis. Moreover, Zhang and Shen [42] studied constraint preconditioners for solving singular saddle point problems, and proved that for solving S2 P2 by the preconditioned GMRES methods with constraint preconditioners, GMRES will determine the least squares solutions at breakdown. Zhang and Wei [39–41] investigated general stationary iterative methods and GMRES methods for solving S2 P2 . In this paper we investigate Uzawa-type methods, which cover the Uzawa method, the preconditioned Uzawa method, and the parameterized Uzawa method for solving singular saddle point problems. Since several Uzawa-type methods are special cases of the abstract inexact Uzawa method [6], we mainly concern ourselves with the semi-convergence of the inexact Uzawa method for solving S2 P2 (1.1). We assume the singular saddle point problem (1.1) satisfying the following assumptions in this paper: • Matrix B is rank-deficient, i.e., rank(B) < m; • Matrix A is positive real, i.e., v T Av > 0, for any nonzero v ∈ Rn ; • (1.1) is consistent, i.e., the solution of (1.1) exists. These assumptions are satisfied when the partial differential equations are discretized by finite difference methods or finite element methods [28,44]. The rest of this paper is organized as follows. In Section 2, we consider how to solve the equivalent linear system of (1.1) and introduce Uzawa-type methods. In Section 3, we make a semi-convergence analysis of Uzawa-type methods for solving singular saddle point problems and discuss the optimal parameters. In Section 4, some numerical examples are given to demonstrate the semi-convergence of Uzawa-type methods. 2. Uzawa-type methods To solve (1.1), we can solve the equivalent linear system instead of the original linear system:
BT 0
A −B
BT 0
A −B
x y
=
f
−g
.
(2.1)
Denote
A≡
,
X≡
x , y
F ≡
f
−g
.
Then (2.1) can be rewritten as
AX = F .
(2.2)
It is well known that the coefficient matrix of (1.1) is indefinite. Whether the coefficient matrix of (1.1) is symmetric or not, the coefficient matrix A of (2.1) is nonsymmetric. Although A loses symmetry, similar to the nonsingular case (cf. [27]), it retains some good properties as follows. Lemma 2.1 ([27]). Assume that A is positive real, then (1) A is range-symmetric, that is, R(A) = R(AT ), or equivalently, N (A) = N (AT ). (2) A is semi-positive real: v T Av ≥ 0, for any v ∈ Rn+m . (3) A is positive semi-stable, i.e., the eigenvalues of A have nonnegative real part.
336
N. Zhang et al. / Journal of Computational and Applied Mathematics 255 (2014) 334–345
For the coefficient matrix A of the linear equations (2.1), we utilize the matrix splitting [35]:
A≡
BT 0
A −B
=
A + Q1 −B
0 Q2
−
Q1 0
−BT Q2
≡ M − P,
(2.3)
where A + Q1 is nonsingular and Q2 is symmetric positive definite, then it yields the inexact Uzawa method [6]: xk+1 = xk + (A + Q1 )−1 (f − Axk − BT yk ), yk+1 = yk + Q2−1 (Bxk+1 − g ),
(2.4)
and the iterative matrix T is:
T = M −1 P = Denote Xk ≡
(A + Q1 )−1 Q1 −1 Q2 B(A + Q1 )−1 Q1
−(A + Q1 )−1 BT . I − Q2−1 B(A + Q1 )−1 BT
(2.5)
k x
yk
, then (2.4) can be rewritten as
Xk+1 = T Xk + M−1 F .
(2.6)
With different choices of Q1 and Q2 , (2.4) covers several Uzawa-type methods as follows: • Uzawa method [3]: Q1 = 0, Q2 = τ1 I (τ > 0), xk+1 = xk + A−1 (f − Axk − BT yk ), yk+1 = yk + τ (Bxk+1 − g ).
(2.7)
• Preconditioned Uzawa method [7]: Q1 = 0, k+1 x = xk + A−1 (f − Axk − BT yk ), k+1 y = yk + Q2−1 (Bxk+1 − g ). • Parameterized Uzawa method, or, the generalized successive overrelaxation (GSOR) method [5,43]: Q1 = Q2 = τ1 Q , where 0 < ω ≤ 1, τ > 0, and Q is symmetric positive definite, k+1 x = (1 − ω)xk + ωA−1 (f − BT yk ), yk+1 = yk + τ Q −1 (Bxk+1 − g ).
(2.8) 1
ω
A − A,
(2.9)
3. Semi-convergence analysis In this section, we investigate the convergence analysis of Uzawa-type methods for solving singular saddle point problems. Let us recall some classical results on stationary [40] or nonstationary [45] iterative methods for solving the singular linear equations. For a matrix B ∈ Rn×n , the smallest nonnegative integer j such that rank(Bj ) = rank(Bj+1 ) is called the index of B, and denoted by j = index(B). Thus index(B) is the size of the largest Jordan block corresponding to the zero eigenvalue of B. Moreover, the Drazin inverse [46] of B is defined as the unique matrix BD which satisfies three matrix equations: BD BBD = BD ,
BD B = BBD ,
Bj+1 BD = Bj .
If j = 1, then the Drazin inverse BD reduces to the group inverse, which denoted by Bg . We now introduce the semi-convergence for solving singular linear systems, which has been studied by several authors (cf. [36,47–49]). Definition 3.1 ([5,50,43]). The iterative method (2.6) is semi-convergent if, for any initial vector X0 ∈ Rm+n , the iteration sequence {Xk } produced by the iterative method (2.6) semi-converges to a solution X∗ of singular linear equations (2.2). Moreover, it holds that X∗ = (I − T )g M −1 F + [I − (I − T )g (I − T )]X0 . For the splitting (2.3) and the associated stationary iteration (2.6), it is well known that if A is nonsingular, then the stationary iteration converges if and only if the spectral radius of the iterative matrix T satisfies ρ(T ) < 1, thus limk→∞ T k = 0. As A is singular, then T has eigenvalue 1 and ρ(T ) cannot be smaller than 1. For the iterative matrix T , we introduce its pseudo-spectral radius ν(T ) [47],
ν(T ) = max{|λ| : λ ∈ σ (T ) \ {1}} where σ (T ) is the set of eigenvalues of matrix T . Now we introduce the necessary and sufficient conditions for the semi-convergence of (2.6). Theorem 3.2 ([50]). The iterative method (2.6) is semi-convergent, if and only if index(I − T ) = 1 and ν(T ) < 1.
(3.1)
N. Zhang et al. / Journal of Computational and Applied Mathematics 255 (2014) 334–345
337
Remark 3.3. If ν(T ) < 1 and index(I − T ) = 1 hold, then T is called semi-convergent ([50, p.198]), while limk→∞ T I − (I − T )g (I − T ). Thus the iteration (2.6) is semi-convergent if and only if T is semi-convergent.
k
=
Remark 3.4. It is known that the magnitude of ν(T ) < 1 determines the convergence rate for {T k }, and the smaller ν(T ) is, the higher asymptotic rate of convergence − ln ν(T ) is. For the iteration (2.6), we require the pseudo-spectral radius of its iterative matrix as small as possible. In Section 4, we discuss the optimal parameters which minimize the pseudo-spectral radius ν(T ). 3.1. The conditions for index(I − T ) = 1 In this section we present the conditions for index(I − T ) = 1. Lemma 3.5 ([40]). Index(I − T ) = 1 holds if and only if, for any 0 ̸= Y ∈ R(A), Y ̸∈ N (AM −1 ). Proof. (⇒) For any 0 ̸= Y ∈ R(A), AM −1 Y ̸= 0, we prove index(I − T ) = 1 by contradiction. Suppose index(I − T ) > 1, then there exists a nonzero X ∈ Rm+n such that M −1 AX ̸= 0 and (M −1 A)2 X = 0. Let Y = AX , then Y ̸= 0, (M−1 A)2 X = M−1 AM−1 Y = 0 ⇒ AM−1 Y = 0, which is in contradiction with 0 ̸= Y ∈ R(A), AM−1 Y ̸= 0. (⇐) If index(I − T ) = 1. We can show that for any 0 ̸= Y ∈ R(A), AM−1 Y ̸= 0. Suppose that there exists a Y = AX ̸= 0 such that AM −1 Y = 0. AX ̸= 0 ⇒ M −1 AX ̸= 0. (M −1 A)2 X = M −1 (AM −1 Y ) = 0, then rank[(M −1 A)2 ] < rank(M −1 A), which contradicts with index(I − T ) = 1. Theorem 3.6. Assume that A is positive real, A + Q1 is nonsingular, and Q2 is symmetric positive definite, then index(I − T ) = 1. Proof. Let X = (ξ T , ηT )T , where ξ ∈ Rn , η ∈ Rm such that
Y = AX =
Aξ + BT η −Bξ
̸= 0.
Notice that
M=
A + Q1 −B
0 Q2
,
M
−1
(A + Q1 )−1 = −1 Q2 B(A + Q1 )−1
0 Q2−1
,
and
AM−1 =
A(A + Q1 )−1 + BT Q2−1 B(A + Q1 )−1 −B(A + Q1 )−1
BT Q2−1 0
,
then
(A + BT Q2−1 B)(A + Q1 )−1 (Aξ + BT η) − BT Q2−1 Bξ AM Y = . −B(A + Q1 )−1 (Aξ + BT η) −1
Denote φ = (A + Q1 )−1 (Aξ + BT η), then
AM−1 Y =
(A + BT Q2−1 B)φ − BT Q2−1 Bξ . −Bφ
To prove index(I − T ) = 1, it is sufficient to prove Y ̸∈ N (AM −1 ) by Lemma 3.5. We fulfill the proof based on the following three cases. 0
Case (a) Aξ + BT η = 0: Y = −Bξ
̸= 0, then Bξ ̸= 0. AM−1 Y =
so is Q2 , then it is easy to see that Y ̸∈ N (AM −1
−1
−BT Q2−1 Bξ 0
, note that Q2 is symmetric positive definite,
).
Case (b) Aξ + B η ̸= 0, Bφ ̸= 0: In this case, it is obvious that Y ̸∈ N (AM −1 ) holds. Case (c) Aξ + BT η ̸= 0, Bφ = 0: T
AM−1 Y =
Aφ − BT Q2−1 Bξ 0
.
To prove Y ̸∈ N (AM −1 ), it is sufficient to prove Aφ − BT Q2−1 Bξ ̸= 0. Next we provide the proof for Aφ − BT Q2−1 Bξ ̸= 0 by contradiction. Assume Aφ − BT Q2−1 Bξ = 0, i.e., φ = A−1 BT Q2−1 Bξ . It follows from Aξ + BT η ̸= 0 that φ ̸= 0, which leads to Bξ ̸∈ N (A−1 BT Q2−1 ).
(3.2)
On the other hand, since Bφ = BA−1 BT Q2 Bξ = 0, we have −1
Bξ ∈ N (BA−1 BT Q2−1 ).
(3.3)
338
N. Zhang et al. / Journal of Computational and Applied Mathematics 255 (2014) 334–345
It follows from (3.2) and (3.3) that rank(A−1 BT Q2−1 ) > rank(BA−1 BT Q2−1 ).
(3.4)
In fact, rank(A−1 BT Q2−1 ) = rank(BT ) = rank(B), and notice that A is positive real, then rank(BA−1 BT Q2−1 ) = rank(BA−1 BT ) = rank(B). The above two equations yield rank(A−1 BT Q2−1 ) = rank(BA−1 BT Q2−1 ), which is in contradiction with (3.4).
3.2. The conditions for ν(T ) < 1
u
Let T
v
=λ
u u n m v , v ̸= 0, where u ∈ C , v ∈ C , that is,
Q1 u − BT v = λ(A + Q1 )u, Q2 v = −λBu + λQ2 v.
(3.5)
λ From (3.5), when λ ̸= 1, it holds that v = λ− Q −1 Bu and 1 2
(Q1 − λA − λQ1 )u =
λ T −1 B Q2 Bu. λ−1
(3.6)
Denote
α=
u∗ Q1 u u∗ u
,
β=
u∗ Au u∗ u
,
γ =
u∗ BT Q2−1 Bu
Multiplying both sides of (3.6) from the left with
λ2 +
u∗ u u∗ u∗ u
.
(3.7)
yields
α γ − 2α − β λ+ = 0. α+β α+β
(3.8)
Lemma 3.7. Suppose that A is positive real, A + Q1 is nonsingular, and Q2 is symmetric positive definite, then u = 0, if and only if λ = 1. Proof. If u = 0, then it follows from (3.5) that (1 −λ)Q2 v = 0, notice that v ̸= 0, or else (uT , v T )T = 0, so λ = 1. Conversely, if λ = 1, then −BT v = Au and Bu = 0. Thus B(−A−1 BT v) = Bu = 0, i.e., (BT v)∗ A−1 (BT v) = 0. Notice that A is positive real, so is A−1 , then BT v = 0, and Au = 0, thus u = 0. Lemma 3.8. Assume that A is positive real, Q1 is symmetric positive semi-definite, and Q2 is symmetric positive definite. If u ∈ N (B), then |λ| < 1. Proof. From (3.5), it holds that (1 − λ)Q2 v = 0, and it follows from Lemma 3.7 that λ ̸= 1, then v = 0. Again from (3.5) α we have (1 − λ)Q1 u = λAu, which yields (α + β)λ = α . Notice that α ≥ 0, and the real part of β is positive, so |λ| = |α+β| < 1. To discuss the conditions for |λ| < 1 with u ̸∈ N (B), we introduce a result about the evaluation of the roots of a quadratic equation as follows. Lemma 3.9 ([51]). Both roots of the quadratic equation x2 −px+q = 0 are less than 1 in modulus if and only if |p−¯pq| < 1−|q|2 , where p¯ is the conjugate number of p. If p and q are real, then the necessary and sufficient conditions for both roots of the quadratic equation x2 − px + q = 0 are less than 1 in modulus is |q| < 1 and |p| < 1 + q [19]. Lemma 3.10. Assume A is positive real, Q1 is symmetric positive semi-definite, and Q2 is symmetric positive definite. If u ̸∈ N (B), then |λ| < 1, if and only if
γ < 2c +
4α c 2 c 2 + d2
,
where β = c + id, or, c and d are the real part and imaginary part of β , respectively, and i =
(3.9)
√ −1.
N. Zhang et al. / Journal of Computational and Applied Mathematics 255 (2014) 334–345
339
Proof. Since u ̸= 0, from Lemma 3.7 it holds that λ ̸= 1, which yields Eq. (3.8)), i.e.,
λ2 +
γ − 2α − β α λ+ = 0. α+β α+β
From Lemma 3.9, we have
γ − 2α − β α 2 γ − 2α − β¯ α . − |λ| < 1 ⇔ <1− α+β α +β α +β α + β¯ By a little algebra, we obtain
|λ| < 1 ⇔ (γ − α)β¯ − βα − β β¯ < |α + β|2 − α 2 , i.e.,
|λ| < 1 ⇔ (α − γ )(c − id) + α(c + id) + c 2 + d2 < |α + c + id|2 − α 2 . Denote s = c 2 + d2 + 2α c, and by some simple algebra operations, we obtain
|λ| < 1 ⇔ |s − γ c + iγ d| < s ⇔ |s − γ c + iγ d|2 < s2 ⇔γ 2 c 2 + γ 2 d2 − 2γ cs < 0. Note that Bu ̸= 0, so γ > 0, then
|λ| < 1 ⇔ γ c 2 + γ d2 − 2cs < 0 ⇔γ (c 2 + d2 ) < 2c (c 2 + d2 + 2α c ) 4α c 2 , ⇔γ < 2c + 2 c + d2 which finishes the proof.
Together with Lemmas 3.7, 3.8 and 3.10, we prove the following result. Theorem 3.11. Assume that A is positive real, Q1 is symmetric positive semi-definite, and Q2 is symmetric positive definite. Then ν(T ) < 1 if and only if, when u ̸∈ N (B), (3.9) holds. Corollary 3.12. Suppose that A and Q2 is symmetric positive definite, and Q1 is symmetric positive semi-definite. Then ν(T ) < 1 if and only if, when u ̸∈ N (B),
γ < 2c + 4α.
(3.10)
3.3. Sufficient conditions for the semi-convergence of Uzawa-type methods The following result gives sufficient conditions for the semi-convergence of the inexact Uzawa method. Theorem 3.13. Assume that A is positive real, Q1 is symmetric positive semi-definite, and Q2 is symmetric positive definite. If
λmax (BT Q2−1 B) < λmin (A + AT ) +
λmin (Q1 )λ2min (A + AT ) , ∥A∥22
(3.11)
then the inexact Uzawa method (2.6) is semi-convergent, where λmax and λmin denote the largest and smallest eigenvalues of the corresponding matrix, respectively. Proof. By Theorems 3.6 and 3.11, it only needs to prove (3.9) when u ̸∈ N (B). Notice that
γ =
u∗ (BT Q2−1 B)u u∗ u
≤ λmax (BT Q2−1 B),
then it follows from (3.11) that
γ < λmin (A + AT ) +
λmin (Q1 )λ2min (A + AT ) . ∥A∥22
It is obvious that
β = c + id =
u∗ Au u∗ u
,
i.e., c 2 + d2 = |β|2 ≤ ∥A∥22 ,
(3.12)
340
N. Zhang et al. / Journal of Computational and Applied Mathematics 255 (2014) 334–345
and c=
1 u∗ (A + AT )u u∗ u
2
≥
1 2
λmin (A + AT ).
Together with the above two inequalities and α =
u∗ Q1 u u∗ u
≥ λmin (Q1 ), we have
4α c 2 λmin (Q1 )λ2min (A + AT ) ≤ 2c + . c 2 + d2 ∥A∥22
λmin (A + AT ) +
4α c 2 , c 2 +d2
From (3.12) and (3.13), it is clear that γ < 2c +
which completes the proof.
(3.13)
Corollary 3.14. Suppose that A and Q2 are symmetric positive definite, Q1 is symmetric positive semi-definite, and
λmax (BT Q2−1 B) < 2λmin (A) + 4λmin (Q1 )
λ2min (A) . λ2max (A)
Then the inexact Uzawa method (2.6) is semi-convergent. The following results give sufficient conditions for the semi-convergence of the Uzawa method. Theorem 3.15. Assume that A is positive real, then the Uzawa method (2.7) is semi-convergent, if 0<τ <
1
. λmax (A + AT )−1 BT B
(3.14)
Proof. By Theorems 3.6 and 3.11, we only need to prove (3.9) when u ̸∈ N (B). Notice that
α = 0,
β=
u∗ Au u∗ u
,
γ =τ
u∗ BT Bu u∗ u
,
then (3.9) ⇔ τ
u∗ BT Bu u∗ u
<
u∗ (A + AT )u
u∗ u ⇔ u (τ B B − (A+ AT ))u < 0 ∗
T
1 1 1 1 ⇔ u∗ (A + AT ) 2 I − τ (A + AT )− 2 BT B(A + AT )− 2 (A + AT ) 2 u > 0. 1
1
The above inequality holds if τ λ[(A + AT )− 2 BT B(A + AT )− 2 ] < 1. In fact, 1
1
λ[(A + AT )− 2 BT B(A + AT )− 2 ] = λ[(A + AT )−1 BT B] ≤ λmax [(A + AT )−1 BT B], then from (3.14) it is easy to see that (3.9) holds.
Corollary 3.16. Assume that A is symmetric positive definite, then the Uzawa method (2.7) is semi-convergent, if 0<τ <
2
λmax (A−1 BT B)
.
(3.15)
The following results give sufficient conditions for the semi-convergence of the preconditioned Uzawa method. Theorem 3.17. Assume that A is positive real, Q2 is symmetric positive definite, then the preconditioned Uzawa method (2.8) is semi-convergent, if
λmax [(A + AT )−1 BT Q2−1 B] < 1. Proof. The proof is similar to that of Theorem 3.15.
(3.16)
Corollary 3.18. Suppose that A and Q2 are symmetric positive definite, then the preconditioned Uzawa method (2.8) is semiconvergent, if
λmax (A−1 BT Q2−1 B) < 2.
(3.17)
Theorem 3.19 gives sufficient conditions for the semi-convergence of the parameterized Uzawa method, which is presented in [43], we give another proof for completeness.
N. Zhang et al. / Journal of Computational and Applied Mathematics 255 (2014) 334–345
341
Theorem 3.19 ([43]). Assume that A and Q are symmetric positive definite, then the parameterized Uzawa method (2.9) is semiconvergent, if 0 < ω < 2,
0<τ <
2(2 − ω)
.
ωλmax (Q −1 BA−1 BT )
(3.18)
Lemma 3.20. Suppose that A and Q are symmetric positive definite, then for the parameterized Uzawa method (2.9), it holds that ν(T ) < 1, if 0 < ω < 2,
γ < 2c + 4α.
(3.19)
Proof. Recall Eq. (3.7), and c is the real part of β ,
α=
1 − ω u∗ Au
ω
u∗ u
,
β=c=
u∗ Au u∗ u
,
γ =τ
u∗ BT Q −1 Bu u∗ u
.
If u = 0, then from Lemma 3.7 we obtain λ = 1.
α If u ̸= 0, and u ∈ N (B), then from the proof of Lemma 3.8 we obtain λ = α+ = 1 − ω. If 0 < ω < 2, then it holds that c |λ| < 1.
If u ̸∈ N (B), taking the same proof of Lemma 3.10, then we obtain |λ| < 1 from γ < 2c + 4α , which finishes the proof.
Proof of Theorem 3.19. By Theorem 3.6 and Lemma 3.20, we only need to prove γ < 2c + 4α when u ̸∈ N (B). Denote 2(2−ω) s = ω , then
γ < 2c + 4α ⇔ u∗ (sA − τ BT Q −1 B)u > 0 1 1 1 1 ⇔ u∗ A 2 (sI − τ A− 2 BT Q −1 BA− 2 )A 2 u > 0. 1
1
1
1
The above inequality holds if τ λmax (A− 2 BT Q −1 BA− 2 ) < s, and λmax (A− 2 BT Q −1 BA− 2 ) = λmax (Q −1 BA−1 BT ), then from (3.18), it is easy to see that γ < 2c + 4α holds, which completes the proof. 3.4. The optimal parameters of Uzawa-type methods Among the Uzawa-type methods we considered, both the Uzawa method (2.7) and the parameterized Uzawa method (2.9) have parameters to be chosen. We now discuss the optimal parameters which minimize the pseudo-spectral radius ν(T ). For the parameterized Uzawa method, Zheng, Bai and Yang [43] gave the optimal parameters ωopt and τopt , as shown by the following theorem. Theorem 3.21 ([43]). Let A and Q be symmetric positive definite. Denote the smallest and the largest nonzero eigenvalues of the matrix Q −1 BA−1 BT by λmin and λmax , respectively. Then for the parameterized Uzawa method (2.9), the optimal parameters ωopt and τopt which minimize ν(T ) are
√ 4 λmin λmax , ωopt = √ √ ( λmin + λmax )2
1
τopt = √ , λmin λmax
(3.20)
and
√ √ λmax − λmin min ν(T ) = √ . √ λmax + λmin
(3.21)
The following theorem presents the optimal parameter τopt for the Uzawa method (2.7) to solve S2 P2 (1.1), and the result is similar to that of the nonsingular SP2 . Theorem 3.22. Assume A is symmetric positive definite. Denote the smallest and the largest nonzero eigenvalues of the matrix A−1 BT B by λmin and λmax , respectively. Then for the Uzawa method (2.7), the optimal parameter τopt which minimizes ν(T ) is
τopt =
2
,
(3.22)
λmax − λmin . λmax + λmin
(3.23)
λmin + λmax
and min ν(T ) =
342
N. Zhang et al. / Journal of Computational and Applied Mathematics 255 (2014) 334–345
Proof. Lemma 3.7 implies u = 0 if and only if λ = 1, which has nothing to do with ν(T ). Assume u ̸= 0 and recall (3.7) and (3.8), we have
α = 0,
β=c=
u∗ Au u∗ u
,
γ =τ
u∗ BT Bu u∗ u
,
and
λ2 +
γ −c c
λ = 0.
Then λ = 0 or λ = 1 − Denote θ =
u∗ BT Bu , u∗ Au
γ
. Notice that γ ̸= 0, or else it yields λ = 1, so Bu ̸= 0.
c
then λ = 1 − τ θ . Since 1
1
1
1
uˆ ∗ (A− 2 BT BA− 2 )ˆu (A− 2 uˆ )∗ BT B(A− 2 uˆ ) = , 0<θ = = 1 1 uˆ ∗ uˆ uˆ ∗ uˆ (A 2 u)∗ (A 2 u) u∗ BT Bu
1
1
1
where uˆ = A 2 u, and notice λ(A− 2 BT BA− 2 ) = λ(A−1 BT B), then λmin ≤ θ ≤ λmax . Now, for the following ‘‘min–max’’ problem min τ
max
λmin ≤θ≤λmax
|λ|,
it is easy to see that the optimal parameter τopt is attained when 1 − τ λmax = −(1 − τ λmin ), that is,
τopt =
2
λmin + λmax
,
and
min ν(T ) = 1 − τopt λmin =
λmax − λmin . λmax + λmin
Remark 3.23. In this subsection we discuss the optimal parameters which minimize the pseudo-spectral radius ν(T ). Notice that if Q is the identity matrix, then λmin and λmax are the same as Theorems 3.21 and 3.22, and it is easy to check that
√ √ λmax − λmin λmax − λmin . < √ √ λ λmax + λmin max + λmin
This means min ν(T ) of Theorem 3.21 is smaller than that of Theorem 3.22, or the parameterized Uzawa method is faster than the Uzawa method. 4. Numerical experiments In this section we compare the performance of several Uzawa-type methods. All the computations are implemented in MATLAB on a PC computer with Genuine Intel (R) CPU 1.83 GHz, and 512 MB memory. Consider the Navier–Stokes equations,
−ν △ u + (u · ▽)u + ▽p = f −▽·u=0
in Ω , in Ω ,
(4.1)
where Ω is an open bounded domain in R2 or R3 , the vector field u represents the velocity in Ω , p denotes pressure, and the scalar ν is the viscosity, which is inversely proportional to the Reynolds number. Linearizing Eq. (4.1) by Picard iteration, we get the Oseen equations,
−ν △ u + (w · ▽)u + ▽p = f −▽·u=0
in Ω , in Ω ,
(4.2)
where w satisfies such that ▽ · w = 0. The test problem is a ‘‘leaky’’ two-dimensional lid-driven cavity problem in the square domain: Ω = (0 < x < 1 : 0 < y < 1). Let u = (u, v)T denote the velocity field and w = (a, b)T the wind. The boundary conditions are u = v = 0 on the three fixed walls (x = 0, y = 0, x = 1), and u = 1, v = 0 on the moving wall (y = 1). We take constant ‘‘wind’’ a = 1, b = 2. To discretize (4.2), we use the ‘‘marker and cell’’ (MAC) finite difference scheme [28,52]. Divide Ω into a uniform n × n grid of cells of width h = 1/n. The discrete velocities and pressures are defined on a staggered grid in which the discrete values of u lie at the centers of the cell boundaries orthogonal to the x-axis, the discrete values of v lie at the centers of the cell boundaries orthogonal to the y-axis, and the discrete pressures lie at the cell centers. We obtain the matrix representation of the Oseen equations (4.2),
AX ≡
A B
BT 0
u p
=
f , 0
(4.3)
N. Zhang et al. / Journal of Computational and Applied Mathematics 255 (2014) 334–345
343
Table 1 CPU time, iteration numbers, and residual with θ = 0.01. Iterative method
CPU (s)
IT
RES
Preconditioned GMRES Uzawa Preconditioned Uzawa Parameterized Uzawa
11.43
109 Divergent 399 Divergent
9.13e-09
15.91
9.88e-09
Table 2 CPU time, iteration numbers, and residual with θ = 0.001. Iterative method
CPU (s)
IT
RES
Preconditioned GMRES Uzawa Preconditioned Uzawa Parameterized Uzawa
9.30 8.63 13.64 61.65
38 72 275 500
6.90e-09 9.42e-09 9.90e-09 3.44e-03
where
A=
F1 0
0 F2
∈ R2n(n−1)×2n(n−1) ,
2 B = (B1 , B2 ) ∈ Rn ×2n(n−1) ,
Fi = ν Ai + Ni ∈ Rn(n−1)×n(n−1) ,
(i = 1, 2), u ∈ R2n(n−1) , p ∈ Rn . For convenience we denote l = 2n(n − 1), m = n2 in this section. The coefficient matrix of linear systems (4.3) satisfies the 2
following properties: A is nonsymmetric and positive real, rank(B) = m−1, rank(A) = l+m−1, N (BT ) = span{em }, em = (1, 1, . . . , 1)T ∈ Rm . To ensure the (1, 1)-block matrix is symmetric positive definite, finally we take the test coefficient matrix as follows,
A¯ =
BT 0
A¯ −B
=
1
(A + AT ) 2 −B
BT 0
.
(4.4)
We compare Uzawa-type methods with the preconditioned GMRES method, and the preconditioner is the block triangular preconditioner [28],
¯
A
Q=
0
BT 1 2 h I
.
(4.5)
ν
For the preconditioned Uzawa method (2.8) and the parameterized Uzawa method (2.9), Q2 and Q are taken as the identity matrix I. Let n = 32, ν = 0.5, then the smallest and the largest nonzero eigenvalues of the matrix A¯ −1 BT B are about 0.0151 and 0.0625, respectively. We take the optimal parameters by Theorems 3.21 and 3.22. For all iterative methods, the initial vector is the zero vector. When the iteration number is k, denote the k-th iterative residual by r k , the stopping criteria is:
∥r k ∥2 ≤ 10−8 . ∥r 0 ∥2 We use the preconditioned restarted GMRES(20) in our numerical tests. The iterative number k of GMRES means k = (i − 1) × 20 + j, where i is the number of restarting, and j is the iteration number of the last restarting, respectively. Notice that it is necessary to compute A¯ −1 for all Uzawa-type methods, and for the preconditioned GMRES, to compute the matrix–vector products with Q−1 , it is still necessary to compute A¯ −1 . Thus, for all iterative methods, we use the MATLAB code cholinc (A¯ , θ ) to approximate A¯ −1 , where θ is the dropping tolerance. With respect to different dropping tolerance θ , we list the CPU time (denoted as CPU), iteration numbers (denoted as IT), and the norm of the residual (denoted as RES) for the above different iterative methods in the tables. Remark 4.1. It can be seen from Tables 1–4, for the Uzawa method, especially for the parameterized Uzawa method, the computation results depend on the exact computation of A¯ −1 , and for the preconditioned GMRES and the preconditioned Uzawa, the computational results do not rely on the exact computation of A¯ −1 . Moreover, the preconditioned GMRES and the preconditioned Uzawa do not require the estimation of the smallest and the largest nonzero eigenvalues of the matrix A¯ −1 BT B. Acknowledgments We would like to express our sincere thanks to our editor Prof. Michael K. Ng and three anonymous reviewers for their valuable suggestions that greatly improved the representation of this paper. Finally, we thank Prof. Eric King-wah Chu for carefully reading our draft and his many valuable suggestions.
344
N. Zhang et al. / Journal of Computational and Applied Mathematics 255 (2014) 334–345 Table 3 CPU time, iteration numbers, and residual with θ = 0.0001. Iterative method
CPU (s)
IT
RES
Preconditioned GMRES Uzawa Preconditioned Uzawa Parameterized Uzawa
8.75 7.67 10.90 61.41
26 22 153 500
1.92e-09 8.96e-09 9.96e-09 6.99e-04
Table 4 CPU time, iteration numbers, and residual with θ = 0.00001. Iterative method
CPU (s)
IT
RES
Preconditioned GMRES Uzawa Preconditioned Uzawa Parameterized Uzawa
8.72 7.64 7.84 61.55
24 16 23 500
5.54e-09 8.89e-09 9.75e-09 1.49e-05
The first author was supported by Zhejiang Provincial Natural Science Foundation of China under grant No. Y1110451 and National Natural Science Foundation of China under grant No. 61002039. The third author was supported by the National Natural Science Foundation of China under grant 11271084, the Doctoral Program of the Ministry of Education under grant 20090071110003, 973 Program Project (No. 2010CB327900) and Shanghai Education Committee. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38]
Z.-Z. Bai, Structured preconditioners for nonsingular matrices of block two-by-two structures, Math. Comp. 75 (2006) 791–815. M. Benzi, G.H. Golub, J. Liesen, Numerical solution of saddle point problems, Acta Numer. 14 (2005) 1–137. K. Arrow, L. Hurwicz, H. Uzawa, Studies in Nonlinear Programming, Stanford University Press, Stanford, CA, 1958. Z.-Z. Bai, Z.-Q. Wang, On parameterized inexact Uzawa methods for generalized saddle point problems, Linear Algebra Appl. 428 (2008) 2900–2932. Z.-Z. Bai, B.N. Parlett, Z.-Q. Wang, On generalized successive overrelaxation methods for augmented linear systems, Numer. Math. 102 (2005) 1–38. 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. H.C. Elman, G.H. Golub, Inexact and preconditioned Uzawa algorithms for saddle point problems, SIAM J. Numer. Anal. 31 (1994) 1645–1661. Q. Hu, J. Zou, An iterative method with variable relaxation parameters for saddle point problems, SIAM J. Matrix Anal. Appl. 23 (2001) 317–338. Q. Hu, J. Zou, Two new variants of nonlinear inexact Uzawa algorithms for saddle point problems, Numer. Math. 93 (2002) 333–359. Y. Lin, Y. Wei, Fast corrected Uzawa methods for solving symmetric saddle point problems, Calcolo 43 (2006) 65–82. J. Lu, Z. Zhang, A modified nonlinear inexact Uzawa algorithm with a variable relaxation parameter for the stabilized saddle point problem, SIAM J. Matrix Anal. Appl. 31 (2010) 1934–1957. W. Zulehner, Analysis of iterative methods for saddle point problems: a unified approach, Math. Comp. 71 (2002) 479–505. 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, 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. 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. Z.-Z. Bai, Optimal parameters in the HSS-like methods for saddle point problems, Numer. Linear Algebra Appl. 16 (2009) 447–479. M.-Q. Jiang, Y. Cao, On local Hermitian and skew-Hermitian splitting iteration methods for generalized saddle point problems, J. Comput. Appl. Math. 231 (2009) 973–982. X. Peng, W. Li, The alternating-direction iterative method for saddle point problems, Appl. Math. Comput. 216 (2010) 1845–1858. D.M. Young, Iterative Solution for Large Linear Systems, Academic Press, New York, 1971. G.H. Golub, X. Wu, J.-Y. Yuan, SOR-like methods for augmented systems, BIT 41 (2001) 71–85. X. Shao, H. Shen, C. Li, The generalized SOR-like method for the augmented systems, Int. J. Inf. Syst. Sci. 2 (1) (2006) 92–98. 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, 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. Z.-Z. Bai, M.K. Ng, On inexact preconditioners for nonsymmetric matrices, SIAM J. Sci. Comput. 26 (2005) 1710–1724. Z.-Z. Bai, M.K. Ng, Z.-Q. Wang, Constraint preconditioners for symmetric indefinite matrices, SIAM J. Matrix Anal. Appl. 31 (2009) 410–433. Z.-Z. Bai, G.H. Golub, C.-K. Li, Convergence properties of preconditioned Hermitian and skew-Hermitian splitting methods for non-Hermitian positive semidefinite matrices, Math. Comp. 76 (2007) 287–298. M. Benzi, G.H. Golub, A preconditioner for generalized saddle point problems, SIAM J. Matrix Anal. Appl. 26 (2004) 20–41. H.C. Elman, Preconditioning for the steady-state Navier–Stokes equations with low viscosity, SIAM J. Sci. Comput. 20 (4) (1999) 1299–1316. C. Keller, N.I.M. Gould, A.J. Wathen, Constraint preconditioning for indefinite linear systems, SIAM J. Matrix Anal. Appl. 21 (2000) 1300–1317. Y. Lin, Y. Wei, A note on constraint preconditioners for nonsymmetric saddle point problems, Numer. Linear Algebra Appl. 14 (2007) 659–664. J.-Y. Pan, M.K. Ng, Z.-Z. Bai, New preconditioners for saddle point problems, Appl. Math. Comput. 172 (2006) 762–771. E.D. Sturier, J. Liesen, Block-diagonal preconditioners for indefinite linear algebraic systems. Part I: theory, SIAM J. Sci. Comput. 26 (5) (2005) 1598–1619. N.I.M. Gould, M.E. Hribar, J. Nocedal, On the solution of equality constrained quadratic programming problems arising in optimization, SIAM J. Sci. Comput. 23 (2001) 1376–1395. V. Sarin, A. Sameh, An efficient iterative method for the generalized Stokes problem, SIAM J. Sci. Comput. 19 (1998) 206–226. F. Chen, Y.-L. Jiang, A generalization of the inexact parameterized Uzawa methods for saddle point problems, Appl. Math. Comput. 206 (2008) 765–771. Z.-Z. Bai, On semi-convergence of Hermitian and skew-Hermitian splitting methods for singular linear systems, Computing 89 (2010) 171–197. Z. Cao, Comparison of performance of iterative methods for singular and nonsingular saddle point linear systems arising from Navier–Stokes equations, Appl. Math. Comput. 174 (2006) 630–642. H. Ma, N. Zhang, A note on block-diagonally preconditioned PIU methods for singular saddle point problems, Int. J. Comput. Math. 88 (2011) 3448–3457.
N. Zhang et al. / Journal of Computational and Applied Mathematics 255 (2014) 334–345
345
[39] N. Zhang, Y. Wei, Solving EP singular linear systems, Int. J. Comput. Math. 81 (2004) 1395–1405. [40] N. Zhang, Y. Wei, On the convergence of general stationary iterative methods for Range–Hermitian singular linear systems, Numer. Linear Algebra Appl. 17 (2010) 139–154. [41] N. Zhang, A note on preconditioned GMRES for solving singular linear equations, BIT 50 (2010) 207–220. [42] N. Zhang, P. Shen, Constraint preconditioners for solving singular saddle point problems, J. Comput. Appl. Math. 238 (2013) 116–125. [43] B. Zheng, Z.-Z. Bai, X. Yang, On semi-convergence of parameterized Uzawa methods for singular saddle point problems, Linear Algebra Appl. 431 (2009) 808–817. [44] X. Wu, B.P.B. Silva, J.Y. Yuan, Conjugate gradient method for rank deficient saddle point problems, Numer. Algorithms 35 (2004) 139–154. [45] X. Shi, Y. Wei, W. Zhang, Convergence of general nonstationary iterative methods for solving singular linear equations, SIAM J. Matrix Anal. Appl. 32 (2011) 72–89. [46] S. Campbell, C. Meyer, Generalized Inverses of Linear Transformations, Pitman, London, 1979, Dover, New York, 1991; SIAM, Philadelphia, 2008. [47] Z. Cao, On the convergence of iterative methods for solving singular linear systems, J. Comput. Appl. Math. 145 (2002) 1–9. [48] W. Li, Y. Liu, X. Peng, The generalized HSS method for solving singular linear systems, J. Comput. Appl. Math. 236 (2012) 2338–2353. [49] Y. Song, L. Wang, Semiconvergence of P-regular splittings for solving singular linear systems, Calcolo 45 (2008) 247–261. [50] A. Berman, R. Plemmons, Nonnegative Matrices in Mathematical Science, Academic Press, New York, 1979, SIAM, 1994. [51] 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. [52] F. Harlow, J. Welch, Numerical calculation of times-dependent viscous incompressible flow of fluid with free surface, Phys. Fluids 8 (1965) 2182–2189.