Volume 2009, Article ID 468965,13pages doi:10.1155/2009/468965
Research Article
New Block Triangular Preconditioners for Saddle Point Linear Systems with
Highly Singular (1,1) Blocks
TingZhu Huang, GuangHui Cheng, and Liang Li
School of Applied Mathematics/Institue of Computational Science, University of Electronic Science and Technology of China, Chengdu, Sichuan 610054, China
Correspondence should be addressed to GuangHui Cheng,[email protected] Received 10 December 2008; Revised 29 March 2009; Accepted 6 August 2009 Recommended by Victoria Vampa
We establish two types of block triangular preconditioners applied to the linear saddle point problems with the singular1,1block. These preconditioners are based on the results presented in the paper of Rees and Greif2007. We study the spectral characteristics of the preconditioners and show that all eigenvalues of the preconditioned matrices are strongly clustered. The choice of the parameter is involved. Furthermore, we give the optimal parameter in practical. Finally, numerical experiments are also reported for illustrating the efficiency of the presented preconditioners.
Copyrightq2009 TingZhu Huang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
Consider the following saddle point linear system:
Aξ≡
G BT B 0
x y
b q
≡f, 1.1
where G ∈ Rn×n is a symmetric and positive semidefinite matrix with nullity dim kernelGp, the matrixB∈Rm×nhas full row rank, vectorsx, b∈Rn,and vectorsy, q∈Rm, and vectorsx, yare unknown. The assumption thatAis nonsingular implies that nullG∩ nullB {0}, which we use in the following analysis. Under these assumptions, the system 1.1has a unique solution. This system is very important and appears in many different applications of scientific computing, such as constrained optimization 1, 2, the finite element method for solving the Navier-Stokes equation3–6, fluid dynamics, constrained least problems and generalized least squares problems 7–10, and the discretized time- harmonic Maxwell equations in mixed form11.
Recently, T. Rees and C. Greif explored a preconditioning technique applied to the problem of solving linear systems arising from primal-dual interior point algorithms and quadratic programming in 12. The preconditioner has the attractive property of improved eigenvalue clustering with increasing ill-conditioned1,1block of the symmetric saddle point systems. To solve the saddle point system1.1, Krylov subspace methods are usually used in modern solution techniques which rely on the ease of sparse matrix-vector products and converges at a rate dependent on the number of distinct eigenvalues of the preconditioned matrix13,14.
The rest of this paper, two types of block triangular preconditioners are established for the saddle point systems with an ill-conditioned1,1block. Our methodology extends the recent work done by Greif and Sch ¨otzau11,15, and Rees and Greif12.
This paper is organized as follows. InSection 2, we will establish new precondtioners and study the spectral analysis of the new preconditioners for the saddle point system. Some numerical examples are given inSection 3. Finally, conclusions are made inSection 4.
2. Preconditioners and Spectrum Analysis
For linear systems, the convergence of an applicable iterative method is determined by the distribution of the eigenvalues of the coefficient matrix. In particular, it is desirable that the number of distinct eigenvalues, or at least the number of clusters, is small, because in this case convergence will be rapid. To be more precise, if there are only a few distinct eigenvalues, then optimal methods like BiCGStab or GMRES will terminatein exact arithmeticafter a small and precisely defined number of steps.
Rees and Greif12established the following preconditioner for the symmetric saddle point system1.1:
M
G BTW−1B tBT
0 W
, 2.1
wheretis a scalar andW is anm×msymmetric positive weight matrix. Similar toM, we introduce the following precondtioners for solving symmetric saddle point systems:
Mt
G BTW−1B 1−tBT
0 tW
, 2.2
wheret /0 is a parameter, and
Mt
⎛
⎝G tBTW−1B tBT
0 1−t
t W
⎞
⎠, 2.3
where 1/t >0.
Theorem 2.1. The matrixM−1t Ahas two distinct eigenvalues which are given by
λ11, λ2−1
t 2.4
with algebraic multiplicitynandp, respectively. The remainingm−peigenvalues satisfy the relation
λ −μ
t μ 1, 2.5
whereμare somem−pgeneralized eigenvalues of the following generalized eigenvalue problem:
BTW−1BxμGx. 2.6
Let{zi}n−mi1 be a basis of the null space ofB, let{ui}pi1be a basis of the null space ofG, and{vi}m−pi1 be a set of linearly independent vectors that complete nullG∪ nullB to a basis of Rn. Then the vectorszTi,0TTi 1, . . . , n−m, the vectorsuTi,1/tW−1BuiTTi 1, . . . , p, and the vectorsvTi,1/tW−1BviTTi 1, . . . , m−p, are linearly independent eigenvectors associated withλ 1, and the vectorsuTi,−W−1BuiTTi 1, . . . , pare linearly independent eigenvectors associated withλ−1/t.
Proof. Suppose thatλis an eigenvalue ofM−1t A, whose eigenvector isx
y
. So, we have
M−1t A x
y
λ x
y
. 2.7
Furthermore, it satisfies the generalized eigenvalue problem G BT
B 0 x
y
λ
G BTW−1B 1−tBT
0 tW
x y
. 2.8
The second block row givesy 1/λtW−1Bx, substituting which into the first block row equation gives
λ−1
λGx
λ 1 t
BTW−1Bx
0. 2.9
By inspection it is straightforward to see that any vectorx∈Rnsatisfies2.9withλ1; thus the latter is an eigenvalue ofM−1t AandxT,1/tW−1BxTTis an eigenvector ofM−1t A. We obtain that the eigenvalueλ1 has algebraic multiplicityn. From the nullity ofGit follows that there areplinearly independent null vectors ofG. For each such null vectorx∈nullG we can obtain
λ−1
t, 2.10
each with algebraic multiplicitypandxT,−W−1BxTTis an eigenvalue ofM−1t A.
Let the vectors {zi}n−mi1 be a basis of the null space of B, and let {ui}pi1 be a basis of the null space of G. Because nullG ∩ nullB {0}, the vectors {zi}n−mi1 and {ui}pi1are linearly independent and together span the subspace nullG∪ nullB. Let the vectors {vi}m−pi1 complete nullG ∪ nullB to a basis of Rn. It follows that the vectors zTi,0TTi 1, . . . , n−m, the vectorsuTi,1/tW−1BuiTTi 1, . . . , p, and the vectors vTi,1/tW−1BviTTi1, . . . , m−p, are linearly independent eigenvectors associated with λ 1, and the vectorsuTi,−W−1BuiTTi 1, . . . , pare linearly independent eigenvectors associated withλ−1/t.
Next, we consider the remainingm−peigenvalues. Supposeλ /1 andλ /−1/t. From 2.9we obtain
BTW−1BxμGx, 2.11
whereμ−tλ/tλ 1, which implies thatλ−μ/tμ 1.
When the parametert−1, we easily obtain the following corollary fromTheorem 2.1.
Corollary 2.2. Lett−1. Then the matrixM−1t Ahas one eigenvalue which is given byλ1 with algebraic multiplicityn p. The remainingm−peigenvalues satisfy the relation
λ μ
μ 1, 2.12
whereμare somem−pgeneralized eigenvalues of the following generalized eigenvalue problem:
BTW−1BxμGx. 2.13
Theorem 2.3. The matrixM−1t Ahas two distinct eigenvalues which are given by
λ11, λ2−1
t 2.14
with algebraic multiplicitynandp, respectively. The remainingm−peigenvalues lie in the interval
0,−1 t
t <0 or
−1 t,0
t >0. 2.15
Proof. According toTheorem 2.1, we know that the matrixM−1t Ahas two distinct eigenvalues which are given by
λ11, λ2−1
t 2.16
with algebraic multiplicitynandp, respectively.
From2.9, we can obtain that the remainingm−peigenvalues satisfy
λ−u
t, 2.17
where u BTW−1Bx, x/Gx, x BTW−1Bx, x ∈ R , in which ·,· is the standard Euclidean inner product, x /∈ nullG and x /∈ nullB. Evidently, we have 0 < u < 1.
The expression2.17gives an explicit formula in terms of the generalized eigenvalues of 2.17and can be used to identify the intervals in which the eigenvalues lie. Furthermore, we can obtain that the remainingm−p eigenvalues lie in the interval0,−1/tt < 0 or
−1/t,0t >0.
When the parametert−1, we easily obtain the following corollary fromTheorem 2.3.
Corollary 2.4. Lett−1. Then the matrixM−1t Ahas one eigenvalue which is given byλ1 with algebraic multiplicityn p. The remainingm−peigenvalues lie in the interval0,1.
Theorem 2.5. The matrixM−1t Ahas two distinct eigenvalues which are given by
λ11, λ2 1
t−1 2.18
with algebraic multiplicitynandp, respectively. The remainingm−peigenvalues satisfy the relation
λ μt
t−1 μt 1, 2.19
whereμare somem−pgeneralized eigenvalues of the following generalized eigenvalue problem:
BTW−1BxμGx. 2.20
Let{zi}n−mi1 be a basis of the null space ofB, and let{ui}pi1be a basis of the null space ofG, and let {vi}m−pi1 be a set of linearly independent vectors that complete nullG∪ nullBto a basis ofRn. Then the vectorszTi,0TTi1, . . . , n−m, the vectorsuTi,t/1−tW−1BuiTTi1, . . . , p, and the vectorsviT,t/1−tW−1BviTTi1, . . . , m−p, are linearly independent eigenvectors associated withλ 1, and the vectorsuTi,−tW−1BuiTTi 1, . . . , pare linearly independent eigenvectors associated withλ1/t−1.
Proof. The proof is similar to the proof ofTheorem 2.1.
When the parametert2, we easily obtain the following corollary fromTheorem 2.5.
Corollary 2.6. Lett2. Then the matrixM−1t Ahas one eigenvalue which is given by
λ1 2.21
with algebraic multiplicityn p. The remainingm−peigenvalues satisfy the relation
λ 2μ
2μ 1, 2.22
whereμare somem−pgeneralized eigenvalues of the following generalized eigenvalue problem:
BTW−1BxμGx. 2.23
Theorem 2.7. The matrixM−1t Ahas two distinct eigenvalues which are given by
λ11, λ2 1
t−1 2.24
with algebraic multiplicitynandp, respectively. The remainingm−peigenvalues lie in the interval
0, 1 t−1
t >1 or 1
t−1,0
t <1. 2.25
Proof. The proof is similar to the proof ofTheorem 2.3.
When the parametert2, we easily obtain the following corollary fromTheorem 2.7.
Corollary 2.8. Lett 2. Then the matrixM−1t Ahas only one eigenvalue which is given byλ 1 with algebraic multiplicityn p. The remainingm−peigenvalues lie in the interval0,1.
Remark 2.9. The above theorems and corollaries illustrate the strong spectral clustering when the1, 1block ofAis singular. A well-known difficulty is the increasing ill-conditioned1, 1block as the solution is approached. Our claim is that the preconditioners perform robust even as the problem becomes more ill-conditioned; in fact the outer iteration count decreases.
On the other hand, solving the augmented1,1block may be more computationally difficult and requires effective approaches such as inexact solvers. InSection 3, we indeed consider inexact solvers in numerical experiments.
Remark 2.10. It is clearly seen from Theorems2.1and2.5and Corollaries2.2and2.6that our preconditioners are suitable for symmetric saddle point systems, from Theorems2.3and2.7 and Corollaries2.4and2.8that our preconditioners are most effective than the preconditioner of12.
0 10 20 30 40 50 60 Drop tolerance0.01
−7
−6
−5
−4
−3
−2
−1 0
a
0 5 10 15 20 25 30 35 40 45
Drop tolerance0.005
−7
−6
−5
−4
−3
−2
−1 0
b
0 5 10 15 20 25 30 35
Drop tolerance0.001
−7
−6
−5
−4
−3
−2
−1 0
t−1 t1 t2
c
0 5 10 15 20 25 30
Drop tolerance0.0001
−7
−6
−5
−4
−3
−2
−1 0
t−1 t1 t2
d
Figure 1: Convergence curve and total numbers of inner GMRES10iterations for differenttwhenh 1/16.
Remark 2.11. Similarly, the nonsymmetric saddle point linear systems can also obtain the above results.
3. Numerical Experiments
All the numerical experiments were performed with MATLAB 7.0. The machine we have used is a PC-IntelR, CoreTM2 CPU T7200 2.0 GHz, 1024 M of RAM. The stopping criterion isrk2/r02 10−6, whererkis the residual vector afterkth iteration. The right-hand side vectorsbandqare taken such that the exact solutionsxandyare both vectors with all components being 1. The initial guess is chosen to be zero vector. We will use preconditioned GMRES10to solve the saddle point linear systems.
0 10 20 30 40 50 60 70 80 Drop tolerance0.01
−7
−6
−5
−4
−3
−2
−1 0 1
a
0 10 20 30 40 50 60
Drop tolerance0.005
−7
−6
−5
−4
−3
−2
−1 0 1
b
0 5 10 15 20 25 30 35 40 45
Drop tolerance0.001
−7
−6
−5
−4
−3
−2
−1 0 1
t−1 t1 t2
c
0 5 10 15 20 25 30
Drop tolerance0.0001
−7
−6
−5
−4
−3
−2
−1 0 1
t−1 t1 t2
d
Figure 2: Convergence curve and total numbers of inner GMRES10iterations for differenttwhenh 1/24.
Our numerical experiments are similar to those in16. We consider the matrices taken from17with notations slightly changed.
We construct the saddle point-type matrix A from reforming a matrix A of the following form:
A
⎛
⎜⎜
⎝
F1 0 BTu 0 F2 BTv Bu Bv 0
⎞
⎟⎟
⎠, 3.1
whereG≡F
1 0 0 F2
is positive real. The matrixAarises from the discretization by the maker and cell finite difference scheme17of a leaky two-dimensional lid-driven cavity problem in a square domain0≤x≤1; 0≤y≤1. Then the matrixBu, Bvis replaced by a random
0 10 20 30 40 50 60 70 80 90 Drop tolerance0.01
−7
−6
−5
−4
−3
−2
−1 0 1
a
0 10 20 30 40 50 60 70 80
Drop tolerance0.005
−7
−6
−5
−4
−3
−2
−1 0 1
b
0 10 20 30 40 50
Drop tolerance0.001
−7
−6
−5
−4
−3
−2
−1 0 1
t−1
t1 t2
c
0 5 10 15 20 25 30 35
Drop tolerance0.0001
−7
−6
−5
−4
−3
−2
−1 0 1
t−1
t1 t2
d
Figure 3: Convergence curve and total numbers of inner GMRES10iterations for differenttwhenh 1/32.
matrixBwith the same sparsity asBu, Bv, replaced byB1 B1 : m,1 :m−3/2Im, such thatB1 is nonsingular. DenoteB2 B1 : m, m 1 : n, then we have B B1, B2with B1∈ Rm,mandB2 ∈ Rm,n−m. Obviously, the resulted saddle point-type matrix
A ≡
G BT B 0
3.2
satisfies rankBT rankB m.
From the matrixAin3.2we construct the following saddle point-type matrix:
A1 ≡
G1 BT B 0
, 3.3
Table 1: Values of n and m, and order ofA1.
h n m Order ofA1
1/16 480 256 736
1/24 1104 576 1680
1/32 1984 1024 3008
Table 2: Number and time of iterations of GMRES10with preconditionersMtandMfor different drop tolerancesτandtwhenh1/16. Results of preconditionerMlie inside .
τ t−1 Time−1 t1 Time1 t2 Time2
0.01 324 0.1875 549 0.3125 655 0.3438
658 0.4219 765 0.4844 1093 0.6563
0.005 322 0.1875 439 0.2656 544 0.2969
547 0.3594 654 0.4063 985 0.6406
0.001 218 0.1719 329 0.2188 330 0.2188
438 0.3594 652 0.4688 985 0.7031
0.0001 216 0.1719 324 0.2188 326 0.2344
436 0.3750 560 0.5156 873 0.7188
whereG1is constructed fromGby making its firstm/4 rows and columns with zero entries.
Note thatG1is semipositive definite and its nullity ism/4.
In our numerical experiments the matrixWin the augmentation block preconditioners is taken asW Im. During implementation of our augmentation block preconditioners, we need the operation G1 BTB−1u for a given vectoru or, equivalently, need to solve the following equation:G1 BTBv ufor which we use an incomplete LU factorization of G1 BTB LU R with drop tolerance τ. Here mn means number of outer inner iterations. Timetrepresents the corresponding computing timein secondswhen taking the parameter ast.
In the following, we summarize the observations from Tables1,2,3,4,5,6, and7and Figures1,2, and3.
iFrom Tables2–4, we can find that our preconditioners are more efficient than those of12both in number of iterations and time of iterations, especially in the case of the optimal parameter.
iiNumber and time of iterations with the preconditionerM−1smaller than those with the preconditionersM1andM2. In fact,M1is a diagonal preconditioner.
iiiNumber and time of iterations with the preconditionerMtare the smallest when t2.
ivNumber of iterations decreases but the computational cost of incomplete LU factorization increases with decreased τ. Therefore, we should not use the preconditioners with smallτin practical.
vThe eigenvalues ofM−1t A1are strongly clustered. Furthermore, the eigenvalues of M−1−1A1are positive.
Table 3: Number and time of iterations of GMRES10with preconditionersMtandMfor different drop tolerancesτandtwhenh1/24. Results of preconditionerMlie inside .
τ t−1 Time−1 t1 Time1 t2 Time2
0.01 328 1.0469 873 2.6563 879 2.8594
988 3.3125 985 3.2188 13125 4.5469
0.005 326 1.0000 659 2.2969 659 2.1719
763 2.4688 872 2.7969 11109 4.2031
0.001 220 0.9219 438 1.5313 543 1.6875
546 1.8438 660 2.4688 1095 3.8750
0.0001 218 0.8594 327 1.2500 328 1.2656
438 1.8750 656 2.6875 987 4.1563 Table 4: Number and time of iterations of GMRES10with preconditionersMtandMfor different drop tolerancesτ,twhenh1/32. Results of preconditionerMlie inside .
τ t−1 Time−1 t1 Time1 t2 Time2
0.01 330 4.0938 983 11.4688 989 12.1563
1093 12.7344 11103 14.1875 13128 17.4219
0.005 328 3.8906 767 9.1563 874 10.1719
876 10.4375 986 11.9688 12118 16.1719
0.001 322 3.1719 544 6.2813 545 6.3438
652 7.3594 765 9.1719 11105 14.7344
0.0001 218 2.8906 330 4.6875 432 5.0313
439 6.0000 658 8.9688 1094 14.5000 Table 5: Number and time of iterations of GMRES10 with preconditioners Mt for different drop tolerancesτandtwhenh1/16.
τ t1/2 Time1/2 t2 Time2 t4 Time4 t8 Time8
0.01 767 0.4844 219 0.1563 323 0.2031 328 0.2188
0.005 653 0.3750 217 0.1563 220 0.1875 324 0.2031
0.001 438 0.3125 214 0.1406 216 0.1719 217 0.1563
0.0001 432 0.3281 213 0.1406 213 0.1406 213 0.1406
Table 6: Number and time of iterations of GMRES10 with preconditioners Mt for different drop tolerancesτandtwhenh1/24.
τ t1/2 Time1/2 t2 Time2 t4 Time4 t8 Time8
0.01 987 3.36 322 0.87 327 1.03 434 1.30
0.005 879 3.02 219 0.88 323 0.90 328 1.14
0.001 652 2.17 216 0.78 219 0.78 322 1.00
0.0001 438 1.77 214 0.75 215 0.73 216 0.83
Table 7: Number and time of iterations of GMRES10 with preconditioners Mt for different drop tolerancesτ,twhenh1/32.
τ t1/2 Time1/2 t2 Time2 t4 Time4 t8 Time8
0.01 11108 15.06 324 3.39 328 3.9 437 5.14
0.005 989 12.31 321 3.02 324 3.39 330 4.17
0.001 658 8.19 217 2.47 220 2.86 324 3.48
0.0001 439 6.06 214 2.28 216 2.55 217 2.69
4. Conclusion
We have proposed two types of block triangular preconditioners applied to the linear saddle point problems with the singular1,1block. The preconditioners have the attractive property of improved eigenvalues clustering with increasing ill-conditioned1,1block. The choice of the parameter is involved. Furthermore, according to Corollaries2.2,2.4,2.6, and2.8, we give the optimal parameter in practice. Numerical experiments are also reported for illustrating the efficiency of the presented preconditioners.
In fact, our methodology can extend the unsymmetrical case; that is, the1,2block and the2,1block of the saddle point linear system are unsymmetrical.
Acknowledgements
Warm thanks are due to the anonymous referees and editor Professor Victoria Vampa who made much useful and detailed suggestions that helped the authors to correct some minor errors and improve the quality of the paper. This research was supported by 973 Program 2008CB317110, NSFC60973015, 10771030, the Chinese Universities Specialized Research Fund for the Doctoral Program of Higher Education 20070614001, and the Project of National Defense Key Lab.9140C6902030906.
References
1 A. Bj ¨orck, “Numerical stability of methods for solving augmented systems,” in Recent Developments˚ in Optimization Theory and Nonlinear Analysis, Y. Censor and S. Reich, Eds., vol. 204 of Contemporary Mathematics, pp. 51–60, American Mathematical Society, Providence, RI, USA, 1997.
2 S. Wright, “Stability of augmented system factorizations in interior-point methods,” SIAM Journal on Matrix Analysis and Applications, vol. 18, no. 1, pp. 191–222, 1997.
3 H. Elman and D. Silvester, “Fast nonsymmetric iterations and preconditioning for Navier-Stokes equations,” SIAM Journal on Scientific Computing, vol. 17, no. 1, pp. 33–46, 1996.
4 H. C. Elman and G. H. Golub, “Inexact and preconditioned Uzawa algorithms for saddle point problems,” SIAM Journal on Numerical Analysis, vol. 31, no. 6, pp. 1645–1661, 1994.
5 H. C. Elman, D. J. Silvester, and A. J. Wathen, “Iterative methods for problems in computational fluid dynamics,” in Iterative Methods in Scientific Computing, R. Chan, T. F. Chan, and G. H. Golub, Eds., pp.
271–327, Springer, Singapore, 1997.
6 B. Fischer, A. Ramage, D. J. Silvester, and A. J. Wathen, “Minimum residual methods for augmented systems,” BIT Numerical Mathematics, vol. 38, no. 3, pp. 527–543, 1998.
7 C. Li, B. Li, and D. J. Evans, “A generalized successive overrelaxation method for least squares problems,” BIT Numerical Mathematics, vol. 38, no. 2, pp. 347–355, 1998.
8 M. F. Murphy, G. H. Golub, and A. J. Wathen, “A note on preconditioning for indefinite linear systems,” Tech. Rep. SCCM-99-03, Stanford University, Stanford, Calif, USA, 1999.
9 S. G. Nash and A. Sofer, “Preconditioning reduced matrices,” SIAM Journal on Matrix Analysis and Applications, vol. 17, no. 1, pp. 47–68, 1996.
10 B. Fischer, A. Ramage, D. J. Silvester, and A. J. Wathen, “Minimum residual methods for augmented systems,” BIT Numerical Mathematics, vol. 38, no. 3, pp. 527–543, 1998.
11 C. Greif and D. Sch ¨otzau, “Preconditioners for the discretized time-harmonic Maxwell equations in mixed form,” Numerical Linear Algebra with Applications, vol. 14, no. 4, pp. 281–297, 2007.
12 T. Rees and C. Greif, “A preconditioner for linear systems arising from interior point optimization methods,” SIAM Journal on Scientific Computing, vol. 29, no. 5, pp. 1992–2007, 2007.
13 J. W. Demmel, Applied Numerical Linear Algebra, SIAM, Philadelphia, Pa, USA, 1997.
14 Y. Saad, Iterative Methods for Sparse Linear Systems, SIAM, Philadelphia, Pa, USA, 2nd edition, 2003.
15 C. Greif and D. Sch ¨otzau, “Preconditioners for saddle point linear systems with highly singular1, 1 blocks,” Electronic Transactions on Numerical Analysis, vol. 22, pp. 114–121, 2006.
16 Z.-H. Cao, “Augmentation block preconditioners for saddle point-type matrices with singular1,1 blocks,” Numerical Linear Algebra with Applications, vol. 15, no. 6, pp. 515–533, 2008.
17 H. C. Elman, “Preconditioning for the steady-state Navier-Stokes equations with low viscosity,” SIAM Journal on Scientific Computing, vol. 20, no. 4, pp. 1299–1316, 1999.
