Convergence of a Generalized USOR Iterative Method for Augmented Systems

Abstract

Zhang and Shang (2010) have presented the Uzawa-SOR (USOR) algorithm to solve augmented systems. In this paper, we establish a generalized Uzawa-SOR (GUSOR) method for solving augmented systems, which is the extension of the USOR method. We prove the convergence of the proposed method under suitable restrictions on the iteration parameters. Lastly, numerical experiments are carried out and experimental results show that our proposed method with appropriate parameters has faster convergence rate than the USOR method.

1. Introduction

We consider the solution of systems of linear equations of the 2-by-2 block structure as follows: 𝐵𝐴𝐵_𝑇0=𝑏𝑞, (1) where 𝐴∈𝑅^𝑚×𝑚 is a symmetric and positive definite matrix and 𝐵∈𝑅^𝑚×𝑛 is a matrix of full column rank, 𝑥,𝑝∈𝑅^𝑚, 𝑦,𝑞∈𝑅^𝑛, and 𝑚≥𝑛. 𝐵^𝑇 denotes the transpose of the matrix 𝐵. We assume that they are of the appropriate dimensions whenever the zero matrix 𝑂 and the identity matrix 𝐼 are used in this paper. The linear systems (1) appears in many different applications of scientific computing, such as the mixed finite element for incompressible flow problems when some form of pressure stabilization is included in the discretization, constrained optimization [1], computational fluid dynamics, and Stokes problems, of constrained least squares problems, and generalized least squares problems [2–5]; see [6, 7] and references therein.

A large variety of methods for solving linear systems of the form (1) can be found in the literature. Yuan and Iusem [3, 5] presented variants of the SOR method and preconditioned conjugate gradient methods. Golub et al. [8] proposed SOR-like algorithms for solving the augmented systems, which was further accelerated and generalized by GSOR method in [9]. Darvishi and Hessari [10] studied the SSOR method. Zhang and Lu [11] studied a GSSOR (generalized SSOR) method. Recently, Zhang and Shang [12] proposed the Uzawa-SOR method and studied its convergence. Bai and Wang [13] established and studied the parameterized inexact Uzawa (PIU) method for solving the corresponding saddle point problems, which was also discussed convergence conditions for matrix splitting iteration methods in [14].

The remainder of the paper is organized as follows. In Section 2, we establish a generalized Uzawa-SOR (GUSOR) method for solving augmented systems and analyze convergence of the corresponding method in Section 3. Numerical results are presented in Section 4. At last, we give some remarks in Section 5.

2. Generalized Uzawa-SOR (GUSOR) Method

For the sake of simplicity, we rewrite system (1) as 𝐴𝐵−𝐵^𝑇0=𝑏−𝑞, (2) where 𝐴∈𝑅^𝑚×𝑚 is a symmetric and positive definite matrix and 𝐵∈𝑅^𝑚×𝑛 is a matrix of full column rank. Let 𝐴 be decomposed as 𝐴=𝐷−𝐿−𝑈 in which 𝐷 is the diagonal of 𝐴, 𝐿 is the strict lower part of 𝐴, and 𝑈 is the strict upper part of 𝐴 with 𝜔 and 𝜏 being nonzero reals.

To construct the generalized USOR method, we consider the following splitting: 𝒜=𝐴𝐵−𝐵^𝑇0=𝐷−𝐿−𝑈𝐵−𝐵^𝑇0=𝒟−ℒ−𝒰, (3) where 𝐵𝒟=𝐷𝑂𝑂𝑄,ℒ=𝐿𝑂_𝑇,,𝛼𝑄𝒰=𝑈−𝐵𝑂𝛽𝑄 (4)𝑄∈ℝ^𝑛×𝑛 is a prescribed a symmetric positive definite matrix and 𝛼+𝛽=1. Let 𝜔 and 𝜏 be two nonzero reals, let 𝐼_𝑚∈𝑅^𝑚×𝑚 and 𝐼_𝑛∈𝑅^𝑛×𝑛 be the m-by-m and the n-by-n identity matrices, respectively, and let Ω be given parameter matrices of the form Ω=𝜔𝐼_𝑚𝑂𝑂𝜏𝐼_𝑛. (5) Then we consider the following generalized SOR iteration scheme for solving the augmented linear system (2): 𝑥^(𝑘+1)𝑦^(𝑘+1)=(𝒟−Ωℒ)⁻¹)]×𝑥[(𝐼−Ω𝒟+Ω𝒰_(𝑘)𝑦^(𝑘)+(𝒟−Ωℒ)⁻¹Ω𝑏,−𝑞 (6) or equivalently, 𝑥^(𝑘+1)𝑦^(𝑘+1)()𝑥=𝐻𝜔,𝜏,𝛼^(𝑘)𝑦^(𝑘)+𝑀(𝜔,𝜏,𝛼)⁻¹𝑏,−𝑞 (7) where𝐻(𝜔,𝜏,𝛼)≡(𝒟−Ωℒ)⁻¹)]=()[(𝐼−Ω𝒟+Ω𝒰𝐷−𝜔𝐿⁻¹𝑂𝜏𝑄1−𝛼𝜏⁻¹𝐵^𝑇(𝐷−𝜔𝐿)⁻¹1𝑄1−𝛼𝜏⁻¹()𝑂()𝑄=()1−𝜔𝐷+𝜔𝑈−𝜔𝐵1−𝛼𝜏𝐷−𝜔𝐿⁻¹((1−𝜔)𝐷+𝜔𝑈)−𝜔(𝐷−𝜔𝐿)⁻¹𝐵𝜏𝑄1−𝛼𝜏⁻¹𝐵^𝑇(𝐷−𝜔𝐿)⁻¹((1−𝜔)𝐷+𝜔𝑈)𝐼−𝜔𝜏𝑄1−𝛼𝜏⁻¹𝐵^𝑇(𝐷−𝜔𝐿)⁻¹𝐵,𝑀(𝜔,𝜏,𝛼)≡Ω⁻¹()=1𝐷−Ω𝐿𝜔(𝐷−𝜔𝐿)𝑂−𝐵^𝑇1−𝛼𝜏𝜏𝑄. (8)(9)

More precisely, we have the following algorithmic description of this GUSOR method.

Generalized USOR Method. Let 𝑄∈ℝ^𝑛×𝑛 be a prescribed symmetric positive definite matrix. Given initial vectors 𝑥⁽⁰⁾∈𝑅^𝑚 and 𝑦⁽⁰⁾∈𝑅^𝑛, and the relaxed parameters 𝛼,𝜔, and 𝜏 with 𝜔,𝜏≠0. For 𝑘=0,1,2,… until the iteration sequence {((𝑥^(𝑘))^𝑇,(𝑦^(𝑘))^𝑇)^𝑇} converges, compute 𝑥^(𝑘+1)=(𝐷−𝜔𝐿)⁻¹×((1−𝜔)𝐷+𝜔𝑈)𝑥^(𝑘)+𝜔𝑏−𝐵𝑦^(𝑘),𝑦^(𝑘+1)=𝑦^(𝑘)+𝜏𝑄1−𝛼𝜏⁻¹𝐵^𝑇𝑥^(𝑘+1)−𝑞. (10)

Remark 1. When the relaxed parameters 𝛼=0, the GUSOR method reduces to the USOR method, so the GUSOR method is the extention of the USOR method.

3. Convergence of the GUSOR Method

In this section, we will analyze a sufficient condition for parameters 𝜔, 𝜏, and 𝛼 in the generalized Uzawa-SOR (GUSOR) method to solve augmented systems (2). We will use the following notations and definitions. For a vector 𝑥, 𝑥^∗ denotes the complex conjugate transpose of the vector 𝑥. 𝜆_min(𝐻) and 𝜆_max(𝐻) denote the minimum and maximum eigenvalues of the Hermitian matrix 𝐻, respectively, and 𝜌(𝐻) denotes the spectral radius of 𝐻. We also assumed that the parameters of 𝛼,𝜔, and 𝜏 used in this paper are positive real numbers.

Note that the iteration matrix of the proposed methods is 𝐻(𝜔,𝜏,𝛼); therefore, the GUSOR method is convergent if and only if the spectral radius of the matrix 𝐻(𝜔,𝜏,𝛼), defined in (8) is less than one; that is, 𝜌(𝐻(𝜔,𝜏,𝛼))<1.

Let 𝜆 be an eigenvalue of 𝐻(𝜔,𝜏,𝛼) and let 𝑥𝑦 be the corresponding eigenvector. Then we have ()𝑂()𝑄𝑥𝑦1−𝜔𝐷+𝜔𝑈−𝜔𝐵1−𝛼𝜏=𝜆𝐷−𝜔𝐿𝑂−𝜏𝐵^𝑇()𝑄𝑥𝑦,1−𝛼𝜏 (11) or equivalently, ()()1−𝜔−𝜆𝐷𝑥+𝜔𝑈𝑥+𝜆𝜔𝐿𝑥=𝜔𝐵𝑦,𝜆−1)(1−𝛼𝜏𝑄𝑦=𝜆𝜏𝐵^𝑇𝑥. (12)

Lemma 2. Let 𝐴 be a symmetric positive definite matrix and 𝐵 a matrix of full column rank. If 𝜆 is an eigenvalue of the iteration matrix 𝐻(𝜔,𝜏,𝛼), then 𝜆≠1.

Proof. If 𝜆=1, since 𝛼 and 𝜏 are positive real numbers, then from (12), we have 𝑥=−𝐴⁻¹𝐵𝑦 and 𝐵^𝑇𝑥=0. It follows that 𝐵^𝑇𝐴⁻¹𝐵𝑦=0 from the above relations, which leads to 𝑦=0 and 𝑥=0. This is a contradiction to the assumption that 𝑥𝑦 is an eigenvector of the iteration matrix 𝐻(𝜔,𝜏,𝛼).

Lemma 3. Let 𝐴 be a symmetric positive definite matrix, and 𝐵 a matrix of full column rank. If 𝜆 is an eigenvalue of the iteration matrix 𝐻(𝜔,𝜏,𝛼) and 𝑥𝑦 is an eigenvector of the iteration matrix 𝐻(𝜔,𝜏,𝛼) corresponding to the eigenvalue 𝜆, then 𝑥≠0. Moreover, if 𝑦=0, then |𝜆|<1.

Proof. The method of proof is exactly the same as in [12], here we omit the proof of Lemma 3.

Lemma 4 (see [13]). Both roots of the complex quadratic equation 𝜆²−𝜙𝜆+𝜓=0 have modulus less than one if and only if |𝜙−𝜙𝜓|+|𝜓|²<1, where 𝜙 denotes the complex conjugate of 𝜙.

Now we are in the position to establish the convergence of the proposed methods. The following theorem presents a sufficient condition for guaranteeing the convergence of the GUSOR method.

Theorem 5. Let 𝐴=𝐷−𝐿−𝑈 be a symmetric positive definite matrix, 𝐵∈𝑅^𝑚×𝑛 a matrix of full column rank, and let 𝑄∈𝑅^𝑛×𝑛 be a symmetric positive definite matrix. Let 𝑆=𝑖(𝐿^𝑇−𝐿)/2. If 𝛼>0, 0<𝜔<2 and 0<𝜏<𝑚₀/(1+𝛼𝑚₀), where 𝑚₀=2(2−𝜔)min_{1≤𝑗≤𝑛}𝑎_𝑗𝑗𝜆²_min(𝐴)𝜔𝜌𝐵𝑄⁻¹𝐵^𝑇4𝜌(𝑆)²+𝜆²_min(𝐴,) (13) then the proposed method is convergent.

Proof. Since 𝜆≠1 from Lemma 2, we have 𝜆≠1. From (12), we obtain 𝑦=𝜆𝜏()𝑄𝜆−1)(1−𝛼𝜏⁻¹𝐵^𝑇()=𝑥,1−𝜔−𝜆𝐷𝑥+𝜔𝑈𝑥+𝜆𝜔𝐿𝑥𝜆𝜔𝜏(𝜆−1)(1−𝛼𝜏)𝐵𝑄⁻¹𝐵^𝑇𝑥. (14)(15) If 𝐵^𝑇𝑥=0, then 𝑦=0 and we have |𝜆|<1 from Lemma 3.
We now assume that 𝐵^𝑇𝑥≠0. For 𝑥≠0, let 𝑥^∗𝐷𝑥𝑥^∗𝑥𝑥=𝛿,^∗𝐿^𝑇𝑥𝑥^∗𝑥𝑥=𝛼+𝛽𝑖,_∗𝐵𝑄⁻¹𝐵^𝑇𝑥𝑥^∗𝑥=𝛾. (16) Notice that 𝑈^𝑇=𝐿, then 𝑥^∗𝑈^𝑇𝑥𝑥^∗𝑥𝑥=𝛼−𝛽𝑖𝛿>0,𝛾>0,𝛿−2𝛼=^∗𝐴𝑥𝑥^∗𝑥>0. (17) Substituting into (15), we obtain ()=1−𝜔−𝜆𝛿+𝜔𝛼−𝛽𝑖+𝜆𝜔𝛼+𝛽𝑖𝜆𝜔𝜏(𝜆−1)(1−𝛼𝜏)𝛾. (18) We notice 𝛿−2𝛼>0 and 0<𝜔<2, that is, −𝛿+𝜔𝛼+𝑖𝜔𝛽≠0, and after some manipulations, we get 𝜆 satisfies the quadratic equation 𝜆²−𝜙𝜆+𝜓=0, where ()((𝜙=𝜔−2𝛿+𝜏/1−𝛼𝜏))𝜔𝛾+2𝜔𝛽𝑖,()−𝛿+𝜔𝛼+𝑖𝜔𝛽𝜑=𝜔−1𝛿+𝜔𝛽𝑖−𝜔𝛼.−𝛿+𝜔𝛼+𝑖𝜔𝛽 (19) Let 𝜒=(−𝛿+𝜔𝛼)²+𝜔²𝛽². By some calculations, one has =1𝜙𝜑−𝜙||𝜒×()1−𝛼𝜏1−𝛼𝜏)(𝜔−2𝛿+𝜏𝜔𝛾²×(𝛿𝜔−2𝜔𝛼)²+2𝜏𝜔²𝛽𝛾^21/2,𝜑²=1𝜒(𝛿(𝜔−1)−𝜔𝛼)²+𝜔²𝛽². (20) From Lemma 4, we know that roots 𝜆 of the complex quadratic equation (18) satisfy |𝜆|<1 if and only if 𝜙−𝜙𝜑+𝜑²<1. (21) By solving (21) for 𝜏, if 0<𝛼𝜏<1, that is, 0<𝜏<1/𝛼, one obtains 𝜏0<<2()1−𝛼𝜏𝛿−2𝛼²(2−𝜔)𝛿4𝛽²+(𝛿−2𝛼)².𝜔𝛾 (22) Letting 2()𝑚=𝛿−2𝛼²(2−𝜔)𝛿4𝛽²+(𝛿−2𝛼)²,𝜔𝛾 (23) we obtain 𝑚0<𝜏<.1+𝛼𝑚 (24) It is obviously that 0<𝜏<1/𝛼. Since 𝑆 is a Hermitian matrix, the eigenvalues of 𝑆 are real, 0≤𝛽²=𝑥^∗𝑆𝑥/𝑥^∗𝑥≤𝜌(𝑆)². For 𝛿−2𝛼=𝑥^∗𝐴𝑥/𝑥^∗𝑥>0, 𝜔>0, 𝛾>0 and 𝑧²/(𝑎+𝑧²) is an increasing function for 𝑧>0, we obtain that |𝜆|<1 when 𝑚0<𝜏<⁰1+𝛼𝑚₀, (25) where 𝑚₀=2(2−𝜔)min_{1≤𝑗≤𝑛}𝑎_𝑗𝑗𝜆²_min(𝐴)𝜔𝜌𝐵𝑄⁻¹𝐵^𝑇4𝜌(𝑆)²+𝜆²_min(𝐴.) (26) Hence, the theorem is proved.

4. Numerical Experiments

In this section, we provide numerical experiments to examine the feasibility and effectiveness of GUSOR method for solving the saddle point problem (1) and compare the results between the GUSOR method and the USOR method provided in [12]. We report the number of iterations (IT), norm of absolution residual vectors (RES), the elapsed CPU time (CPU), and the spectral radius of corresponding iterative matrix denoted by 𝜌. Here, RES is defined as RES=𝑏−𝐴𝑥^(𝑘)−𝐵^𝑇𝑦^(𝑘)2₂+𝑞−𝐵^𝑇𝑥^(𝑘)2₂‖𝑏‖²₂+𝑞²₂, (27) with (𝑥^𝑇_𝑘,𝑦^𝑇_𝑘)^𝑇 being the final approximate solution, where ‖⋅‖ refers to 𝐿₂-norm. We choose the right-hand vector (𝑏^𝑇,𝑞^𝑇)^𝑇∈𝑅^𝑚+𝑛 such that the exact solution of the augmented linear system (1) is (𝑥^𝑇_∗,𝑦^𝑇_∗)^𝑇=(1,1,…,1)^𝑇∈𝑅^𝑚+𝑛. All numerical experiments are carried out on a PC equipped with Intel Core i3 2.3 GHz CPU and 2.00 GB RAM memory Using MATLAB R2010a.

Example 6 (see [9]). Let the augmented system (1) in which 𝐴&=𝐼⊗𝑇+𝑇⊗𝐼𝑂𝑂𝐼⊗𝑇+𝑇⊗𝐼∈𝑅^{2𝑝2×2𝑝2},𝐵=𝐼⊗𝐹𝐹⊗𝐼∈𝑅^{2𝑝2×2𝑝2},1𝑇=ℎ²tridiag(−1,2,−1)∈𝑅^𝑝×𝑝,1𝐹=ℎtridiag(−1,1,0)∈𝑅^𝑝×𝑝, (28) with ⊗ is the Kronecker product symbol and ℎ=1/(𝑝+1) and 𝑆=tridiag(𝑎,𝑏,𝑐) is a tridiagonal matrix with 𝑆_𝑖,𝑖=𝑏, 𝑆_{𝑖−1,𝑖}=𝑎, 𝑆_𝑖,𝑖+1=𝑐 for appropriate 𝑖.

For this example, 𝑚=2𝑝² and 𝑛=𝑝². Hence, the total number of variables is 𝑚+𝑛=3𝑝². We choose the matrix 𝑄 as an approximation to the matrix 𝐵^𝑇𝐴⁻¹𝐵, according to three cases listed in Table 1.

Table 1 

Choices of the matrix 𝑄.

In our experiments, all runs with respect to both USOR method and GUSOR method are started from the initial vector which is set to the zero vector and terminate if the current iteration satisfies ERR<10⁻⁶. Here, ERR is defined as ERR=𝑥_𝑘−𝑥_∗²+𝑦_𝑘−𝑦_∗²𝑥₀−𝑥_∗²+𝑦₀−𝑦_∗². (29)

In Tables 2, 3, and 4, we list the values of (𝜔,𝜏) which are same as in [12], IT, RES, CPU, and the spectral radii of corresponding iterative matrices for various problem sizes (𝑚,𝑛), respectively. They clearly show that the GUSOR method is more effective than the USOR method on convergence rate, computing speed, and the spectral radii of corresponding iterative matrices. IT and CPU of our proposed method are nearly half of the USOR if 𝑝 is smaller. However, the relaxed parameters 𝛼 of GUSOR method are not optimal values and only lie in the convergence region of the method. The determination of optimum values of the parameters needs further study.

Table 2 

Numerical results for Example 6 in Case I with 𝛼=0.4.

Table 3 

Numerical results for Example 6 in Case II with 𝛼=0.6.

Table 4 

Numerical results for Example 4.1 in Case III with 𝛼=0.4.

Remark 7. When 𝛼=0, in this case, the proposed method is the one in [12]. Through experiment results, we find the optimal relaxation of 𝛼 seems to be about 0.5.

5. Conclusions

In this paper, we propose the GUSOR method for the solution of the saddle point problems and analyze the convergence of GUSOR method. When chosen the relaxed parameters 𝛼, the spectral radii of the iteration matrices, IT and CPU with the proposed method are smaller than those in [12], which is shown through numerical experiments. Particularly, one may discuss how to select the set of optimal parameters for accelerating the convergence of the considered method effectively. The optimal choice of this set of parameters is valuably studied which is our future work.

Acknowledgments

The authors would like to thank the anonymous referees for their helpful comments and advice, which greatly improved the paper. The study was financially supported by the National Natural Science Foundation (nos. 11161041 and 71301111), Chinese Universities Specialized Research Fund for the Doctoral Program (20110185110020).