New Block Triangular Preconditioners for Saddle Point Linear Systems with Highly Singular (1,1) Blocks

Abstract

We establish two types of block triangular preconditioners applied to the linear saddle point problems with the singular (1,1) block. These preconditioners are based on the results presented in the paper of Rees and Greif (2007). 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.

1. Introduction

Consider the following saddle point linear system:

𝒜𝜉≡𝐺𝐵^𝑇𝑥𝑦=𝑏𝐵0𝑞≡𝑓, (1.1) where 𝐺∈ℝ^𝑛×𝑛 is a symmetric and positive semidefinite matrix with nullity (= dim(kernel(𝐺))) 𝑝, the matrix 𝐵∈ℝ^𝑚×𝑛 has full row rank, vectors 𝑥,𝑏∈ℝ^𝑛, and vectors 𝑦,𝑞∈ℝ^𝑚, and vectors 𝑥,𝑦 are unknown. The assumption that 𝒜 is nonsingular implies that null(𝐺)∩null(𝐵)={0}, which we use in the following analysis. Under these assumptions, the system (1.1) has 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 equation [3–6], fluid dynamics, constrained least problems and generalized least squares problems [7–10], and the discretized time-harmonic Maxwell equations in mixed form [11].

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-conditioned (1,1) block of the symmetric saddle point systems. To solve the saddle point system (1.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 matrix [13, 14].

The rest of this paper, two types of block triangular preconditioners are established for the saddle point systems with an ill-conditioned (1,1) block. Our methodology extends the recent work done by Greif and Schötzau [11, 15], and Rees and Greif [12].

This paper is organized as follows. In Section 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 in Section 3. Finally, conclusions are made in Section 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 terminate (in exact arithmetic) after a small and precisely defined number of steps.

Rees and Greif [12] established the following preconditioner for the symmetric saddle point system (1.1):

ℳ=𝐺+𝐵^𝑇𝑊⁻¹𝐵𝑡𝐵^𝑇0𝑊, (2.1) where 𝑡 is a scalar and 𝑊 is an 𝑚×𝑚 symmetric positive weight matrix. Similar to ℳ, we introduce the following precondtioners for solving symmetric saddle point systems:

ℳ_𝑡=𝐺+𝐵^𝑇𝑊⁻¹𝐵(1−𝑡)𝐵^𝑇0𝑡𝑊, (2.2) where 𝑡≠0 is a parameter, and

ℳ_𝑡=𝐺+𝑡𝐵^𝑇𝑊⁻¹𝐵𝑡𝐵^𝑇01−𝑡𝑡𝑊, (2.3) where 1≠𝑡>0.

Theorem 2.1. The matrix ℳ_𝑡⁻¹𝒜 has two distinct eigenvalues which are given by 𝜆₁=1,𝜆₂1=−𝑡 (2.4) with algebraic multiplicity 𝑛 and 𝑝, respectively. The remaining 𝑚−𝑝 eigenvalues satisfy the relation 𝜆=−𝜇𝑡(𝜇+1), (2.5) where 𝜇 are some 𝑚−𝑝 generalized eigenvalues of the following generalized eigenvalue problem: 𝐵^𝑇𝑊⁻¹𝐵𝑥=𝜇𝐺𝑥. (2.6) Let {𝑧_𝑖}^𝑛−𝑚_𝑖=1 be a basis of the null space of 𝐵, let {𝑢_𝑖}^𝑝_𝑖=1 be a basis of the null space of 𝐺, and {𝑣_𝑖}^𝑚−𝑝_𝑖=1 be a set of linearly independent vectors that complete null(𝐺)∪null(𝐵) to a basis of ℝ^𝑛. Then the vectors [𝑧^𝑇_𝑖,0^𝑇]^𝑇(𝑖=1,…,𝑛−𝑚), the vectors [𝑢^𝑇_𝑖,(1/𝑡)(𝑊⁻¹𝐵𝑢_𝑖)^𝑇]^𝑇(𝑖=1,…,𝑝), and the vectors [𝑣^𝑇_𝑖,(1/𝑡)(𝑊⁻¹𝐵𝑣_𝑖)^𝑇]^𝑇(𝑖=1,…,𝑚−𝑝), are linearly independent eigenvectors associated with 𝜆=1, and the vectors [𝑢^𝑇_𝑖,−(𝑊⁻¹𝐵𝑢_𝑖)^𝑇]^𝑇(𝑖=1,…,𝑝) are linearly independent eigenvectors associated with 𝜆=−1/𝑡.

Proof. Suppose that 𝜆 is an eigenvalue of ℳ_𝑡⁻¹𝒜, whose eigenvector is (^𝑥_𝑦). So, we have ℳ_𝑡⁻¹𝒜𝑥𝑦𝑥𝑦=𝜆. (2.7) Furthermore, it satisfies the generalized eigenvalue problem 𝐺𝐵^𝑇𝑥𝐵0𝑦=𝜆𝐺+𝐵^𝑇𝑊⁻¹𝐵(1−𝑡)𝐵^𝑇𝑥0𝑡𝑊𝑦. (2.8) The second block row gives 𝑦=(1/𝜆𝑡)𝑊⁻¹𝐵𝑥, substituting which into the first block row equation gives 1(𝜆−1)[𝜆𝐺𝑥+(𝜆+𝑡)𝐵^𝑇𝑊⁻¹𝐵𝑥]=0. (2.9) By inspection it is straightforward to see that any vector 𝑥∈ℝ^𝑛 satisfies (2.9) with 𝜆=1; thus the latter is an eigenvalue of ℳ_𝑡⁻¹𝒜 and [𝑥^𝑇,(1/𝑡)(𝑊⁻¹𝐵𝑥)^𝑇]^𝑇 is an eigenvector of ℳ_𝑡⁻¹𝒜. We obtain that the eigenvalue 𝜆=1 has algebraic multiplicity 𝑛. From the nullity of 𝐺 it follows that there are 𝑝 linearly independent null vectors of 𝐺. For each such null vector 𝑥∈null(𝐺) we can obtain 1𝜆=−𝑡, (2.10) each with algebraic multiplicity 𝑝 and [𝑥^𝑇,−(𝑊⁻¹𝐵𝑥)^𝑇]^𝑇 is an eigenvalue of ℳ_𝑡⁻¹𝒜.
Let the vectors {𝑧_𝑖}^𝑛−𝑚_𝑖=1 be a basis of the null space of 𝐵, and let {𝑢_𝑖}^𝑝_𝑖=1 be a basis of the null space of 𝐺. Because null(𝐺)∩null(𝐵)={0}, the vectors {𝑧_𝑖}^𝑛−𝑚_𝑖=1 and {𝑢_𝑖}^𝑝_𝑖=1 are linearly independent and together span the subspace null(𝐺)∪null(𝐵). Let the vectors {𝑣_𝑖}^𝑚−𝑝_𝑖=1 complete null(𝐺)∪null(𝐵) to a basis of ℝ^𝑛. It follows that the vectors [𝑧^𝑇_𝑖,0^𝑇]^𝑇(𝑖=1,…,𝑛−𝑚), the vectors [𝑢^𝑇_𝑖,(1/𝑡)(𝑊⁻¹𝐵𝑢_𝑖)^𝑇]^𝑇(𝑖=1,…,𝑝), and the vectors [𝑣^𝑇_𝑖,(1/𝑡)(𝑊⁻¹𝐵𝑣_𝑖)^𝑇]^𝑇(𝑖=1,…,𝑚−𝑝), are linearly independent eigenvectors associated with 𝜆=1, and the vectors [𝑢^𝑇_𝑖,−(𝑊⁻¹𝐵𝑢_𝑖)^𝑇]^𝑇(𝑖=1,…,𝑝) are linearly independent eigenvectors associated with 𝜆=−1/𝑡.
Next, we consider the remaining 𝑚−𝑝 eigenvalues. Suppose 𝜆≠1 and 𝜆≠−1/𝑡. From (2.9) we obtain
𝐵^𝑇𝑊⁻¹𝐵𝑥=𝜇𝐺𝑥, (2.11) where 𝜇=−𝑡𝜆/(𝑡𝜆+1), which implies that 𝜆=−𝜇/(𝑡(𝜇+1)).

When the parameter 𝑡=−1, we easily obtain the following corollary from Theorem 2.1.

Corollary 2.2. Let 𝑡=−1. Then the matrix ℳ_𝑡⁻¹𝒜 has one eigenvalue which is given by 𝜆=1 with algebraic multiplicity 𝑛+𝑝. The remaining 𝑚−𝑝 eigenvalues satisfy the relation 𝜇𝜆=𝜇+1, (2.12) where 𝜇 are some 𝑚−𝑝 generalized eigenvalues of the following generalized eigenvalue problem: 𝐵^𝑇𝑊⁻¹𝐵𝑥=𝜇𝐺𝑥. (2.13)

Theorem 2.3. The matrix ℳ_𝑡⁻¹𝒜 has two distinct eigenvalues which are given by 𝜆₁=1,𝜆₂1=−𝑡 (2.14) with algebraic multiplicity 𝑛 and 𝑝, respectively. The remaining 𝑚−𝑝 eigenvalues lie in the interval 1(0,−𝑡1)(𝑡<0)or(−𝑡,0)(𝑡>0). (2.15)

Proof. According to Theorem 2.1, we know that the matrix ℳ_𝑡⁻¹𝒜 has two distinct eigenvalues which are given by 𝜆₁=1,𝜆₂1=−𝑡 (2.16) with algebraic multiplicity 𝑛 and 𝑝, respectively.
From (2.9), we can obtain that the remaining 𝑚−𝑝 eigenvalues satisfy
𝑢𝜆=−𝑡, (2.17) where 𝑢=⟨𝐵^𝑇𝑊⁻¹𝐵𝑥,𝑥⟩/⟨𝐺𝑥,𝑥⟩+⟨𝐵^𝑇𝑊⁻¹𝐵𝑥,𝑥⟩∈ℝ⁺, in which ⟨⋅,⋅⟩ is the standard Euclidean inner product, 𝑥∉null(𝐺) and 𝑥∉null(𝐵). Evidently, we have 0<𝑢<1. The expression (2.17) gives an explicit formula in terms of the generalized eigenvalues of (2.17) and can be used to identify the intervals in which the eigenvalues lie. Furthermore, we can obtain that the remaining 𝑚−𝑝 eigenvalues lie in the interval (0,−1/𝑡)(𝑡<0) or (−1/𝑡,0)(𝑡>0).

When the parameter 𝑡=−1, we easily obtain the following corollary from Theorem 2.3.

Corollary 2.4. Let 𝑡=−1. Then the matrix ℳ_𝑡⁻¹𝒜 has one eigenvalue which is given by 𝜆=1 with algebraic multiplicity 𝑛+𝑝. The remaining 𝑚−𝑝 eigenvalues lie in the interval (0,1).

Theorem 2.5. The matrix ℳ_𝑡⁻¹𝒜 has two distinct eigenvalues which are given by 𝜆₁=1,𝜆₂=1𝑡−1 (2.18) with algebraic multiplicity 𝑛 and 𝑝, respectively. The remaining 𝑚−𝑝 eigenvalues satisfy the relation 𝜆=𝜇𝑡(𝑡−1)(𝜇𝑡+1), (2.19) where 𝜇 are some 𝑚−𝑝 generalized eigenvalues of the following generalized eigenvalue problem: 𝐵^𝑇𝑊⁻¹𝐵𝑥=𝜇𝐺𝑥. (2.20) Let {𝑧_𝑖}^𝑛−𝑚_𝑖=1 be a basis of the null space of 𝐵, and let {𝑢_𝑖}^𝑝_𝑖=1 be a basis of the null space of 𝐺, and let {𝑣_𝑖}^𝑚−𝑝_𝑖=1 be a set of linearly independent vectors that complete null(𝐺)∪null(𝐵) to a basis of ℝ^𝑛. Then the vectors [𝑧^𝑇_𝑖,0^𝑇]^𝑇(𝑖=1,…,𝑛−𝑚), the vectors [𝑢^𝑇_𝑖,(𝑡/(1−𝑡))(𝑊⁻¹𝐵𝑢_𝑖)^𝑇]^𝑇(𝑖=1,…,𝑝), and the vectors [𝑣^𝑇_𝑖,(𝑡/(1−𝑡))(𝑊⁻¹𝐵𝑣_𝑖)^𝑇]^𝑇(𝑖=1,…,𝑚−𝑝), are linearly independent eigenvectors associated with 𝜆=1, and the vectors [𝑢^𝑇_𝑖,−𝑡(𝑊⁻¹𝐵𝑢_𝑖)^𝑇]^𝑇(𝑖=1,…,𝑝) are linearly independent eigenvectors associated with 𝜆=1/(𝑡−1).

Proof. The proof is similar to the proof of Theorem 2.1.

When the parameter 𝑡=2, we easily obtain the following corollary from Theorem 2.5.

Corollary 2.6. Let 𝑡=2. Then the matrix ℳ_𝑡⁻¹𝒜 has one eigenvalue which is given by 𝜆=1 (2.21) with algebraic multiplicity 𝑛+𝑝. The remaining 𝑚−𝑝 eigenvalues satisfy the relation 𝜆=2𝜇2𝜇+1, (2.22) where 𝜇 are some 𝑚−𝑝 generalized eigenvalues of the following generalized eigenvalue problem: 𝐵^𝑇𝑊⁻¹𝐵𝑥=𝜇𝐺𝑥. (2.23)

Theorem 2.7. The matrix ℳ_𝑡⁻¹𝒜 has two distinct eigenvalues which are given by 𝜆₁=1,𝜆₂=1𝑡−1 (2.24) with algebraic multiplicity 𝑛 and 𝑝, respectively. The remaining 𝑚−𝑝 eigenvalues lie in the interval 1(0,1𝑡−1)(𝑡>1)or(𝑡−1,0)(𝑡<1). (2.25)

Proof. The proof is similar to the proof of Theorem 2.3.

When the parameter 𝑡=2, we easily obtain the following corollary from Theorem 2.7.

Corollary 2.8. Let 𝑡=2. Then the matrix ℳ_𝑡⁻¹𝒜 has only one eigenvalue which is given by 𝜆=1 with algebraic multiplicity 𝑛+𝑝. The remaining 𝑚−𝑝 eigenvalues lie in the interval (0,1).

Remark 2.9. The above theorems and corollaries illustrate the strong spectral clustering when the (1, 1) block of 𝒜 is singular. A well-known difficulty is the increasing ill-conditioned (1, 1) block 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 augmented (1,1) block may be more computationally difficult and requires effective approaches such as inexact solvers. In Section 3, we indeed consider inexact solvers in numerical experiments.

Remark 2.10. It is clearly seen from Theorems 2.1 and 2.5 and Corollaries 2.2 and 2.6 that our preconditioners are suitable for symmetric saddle point systems, from Theorems 2.3 and 2.7 and Corollaries 2.4 and 2.8 that our preconditioners are most effective than the preconditioner of [12].

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-Intel(R), Core(TM)2 CPU T7200 2.0 GHz, 1024 M of RAM. The stopping criterion is ‖𝑟^(𝑘)‖₂/‖𝑟⁽⁰⁾‖₂⩽10⁻⁶, where 𝑟^(𝑘) is the residual vector after 𝑘th iteration. The right-hand side vectors 𝑏 and 𝑞 are taken such that the exact solutions 𝑥 and 𝑦 are both vectors with all components being 1. The initial guess is chosen to be zero vector. We will use preconditioned GMRES(10) to solve the saddle point linear systems.

Our numerical experiments are similar to those in [16]. We consider the matrices taken from [17] with notations slightly changed.

We construct the saddle point-type matrix 𝒜 from reforming a matrix 𝒜 of the following form:

𝐹𝒜=¹0𝐵^𝑇_𝑢0𝐹₂𝐵^𝑇_𝑣𝐵_𝑢𝐵_𝑣0, (3.1) where 𝐺≡(^𝐹10_0𝐹2) is positive real. The matrix 𝒜 arises from the discretization by the maker and cell finite difference scheme [17] of a leaky two-dimensional lid-driven cavity problem in a square domain (0≤𝑥≤1;0≤𝑦≤1). Then the matrix [𝐵_𝑢,𝐵_𝑣] is replaced by a random matrix 𝐵 with the same sparsity as [𝐵_𝑢,𝐵_𝑣], replaced by 𝐵₁=𝐵(1∶𝑚,1∶𝑚)−(3/2)𝐼_𝑚, such that 𝐵₁ is nonsingular. Denote 𝐵₂=𝐵(1∶𝑚,𝑚+1∶𝑛), then we have 𝐵=[𝐵₁,𝐵₂] with 𝐵₁∈ℛ^𝑚,𝑚 and 𝐵₂∈ℛ^{𝑚,𝑛−𝑚}. Obviously, the resulted saddle point-type matrix

𝒜≡𝐺𝐵^𝑇𝐵0 (3.2) satisfies rank(𝐵^𝑇)=rank(𝐵)=𝑚.

From the matrix 𝒜 in (3.2) we construct the following saddle point-type matrix:

𝒜₁≡𝐺₁𝐵^𝑇𝐵0, (3.3) where 𝐺₁ is constructed from 𝐺 by making its first 𝑚/4 rows and columns with zero entries. Note that 𝐺₁ is semipositive definite and its nullity is 𝑚/4.

In our numerical experiments the matrix 𝑊 in the augmentation block preconditioners is taken as 𝑊=𝐼_𝑚. During implementation of our augmentation block preconditioners, we need the operation (𝐺₁+𝐵^𝑇𝐵)⁻¹𝑢 for a given vector 𝑢 or, equivalently, need to solve the following equation: (𝐺₁+𝐵^𝑇𝐵)𝑣=𝑢 for which we use an incomplete LU factorization of (𝐺₁+𝐵^𝑇𝐵)=𝐿𝑈+𝑅 with drop tolerance 𝜏. Here 𝑚(𝑛) means number of outer (inner) iterations. Time(𝑡) represents the corresponding computing time (in seconds) when taking the parameter as 𝑡.

In the following, we summarize the observations from Tables 1, 2, 3, 4, 5, 6, and 7 and Figures 1, 2, and 3.

Table 1 

Values of n and m, and order of 𝒜₁.

Table 2 

Number and time of iterations of GMRES(10) with preconditioners ℳ_𝑡 and ℳ for different drop tolerances 𝜏 and 𝑡 when ℎ=1/16. Results of preconditioner ℳ lie inside [].

Table 3 

Number and time of iterations of GMRES(10) with preconditioners ℳ_𝑡 and ℳ for different drop tolerances 𝜏 and 𝑡 when ℎ=1/24. Results of preconditioner ℳ lie inside [].

Table 4 

Number and time of iterations of GMRES(10) with preconditioners ℳ_𝑡 and ℳ for different drop tolerances 𝜏 , 𝑡 when ℎ=1/32. Results of preconditioner ℳ lie inside [].

Table 5 

Number and time of iterations of GMRES(10) with preconditioners ℳ_𝐭 for different drop tolerances 𝜏 and 𝑡 when ℎ=1/16.

Table 6 

Number and time of iterations of GMRES(10) with preconditioners ℳ_𝐭 for different drop tolerances 𝜏 and 𝐭 when ℎ=1/24.

Table 7 

Number and time of iterations of GMRES(10) with preconditioners ℳ_𝐭 for different drop tolerances 𝜏, 𝑡 when ℎ=1/32.

(a)

(b)

(c)

(d)

Figure 1 

Convergence curve and total numbers of inner GMRES(10) iterations for different 𝑡 when ℎ=1/16.

(a)

(b)

(c)

(d)

Figure 2 

Convergence curve and total numbers of inner GMRES(10) iterations for different 𝑡 when ℎ=1/24.

(a)

(b)

(c)

(d)

Figure 3 

Convergence curve and total numbers of inner GMRES(10) iterations for different 𝑡 when ℎ=1/32.

(i)From Tables 2–4, we can find that our preconditioners are more efficient than those of [12] both in number of iterations and time of iterations, especially in the case of the optimal parameter.(ii)Number and time of iterations with the preconditioner ℳ₋₁ smaller than those with the preconditioners ℳ₁ and ℳ₂. In fact, ℳ₁ is a diagonal preconditioner.(iii)Number and time of iterations with the preconditioner ℳ_𝑡 are the smallest when 𝑡=2.(iv)Number 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.(v)The eigenvalues of ℳ_𝑡⁻¹𝒜₁ are strongly clustered. Furthermore, the eigenvalues of ℳ⁻¹₋₁𝒜₁ are positive.

4. Conclusion

We have proposed two types of block triangular preconditioners applied to the linear saddle point problems with the singular (1,1) block. The preconditioners have the attractive property of improved eigenvalues clustering with increasing ill-conditioned (1,1) block. The choice of the parameter is involved. Furthermore, according to Corollaries 2.2, 2.4, 2.6, and 2.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, the (1,2) block and the (2,1) block of the saddle point linear system are unsymmetrical.

Acknowledgments

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), NSFC (60973015, 10771030), the Chinese Universities Specialized Research Fund for the Doctoral Program of Higher Education (20070614001), and the Project of National Defense Key Lab. (9140C6902030906).