New Preconditioners for Nonsymmetric Saddle Point Systems with Singular (1,1) Block

Abstract

We investigate the solution of large linear systems of saddle point type with singular (1,1) block by preconditioned iterative methods and consider two parameterized block triangular preconditioners used with Krylov subspace methods which have the attractive property of improved eigenvalue clustering with increased ill-conditioning of the (1,1) block of the saddle point matrix, including the choice of the parameter. Meanwhile, we analyze the spectral characteristics of two preconditioners and give the optimal parameter in practice. Numerical experiments that validate the analysis are presented.

1. Introduction

We study preconditioners for general nonsingular linear systems of the type 𝒜𝑢=𝐴𝐵^𝑇𝑥𝑦=𝑓𝑔,𝐶0𝐴∈𝑅^𝑛×𝑛,𝐵,𝐶∈𝑅^𝑚×𝑛,𝑛≥𝑚. (1) Such systems arise in a large number of applications, for example, the (linearized) Navier-Stokes equations and other physical problems with conservation laws as well as constrained optimization problems [1–5].

As such systems are typically large and sparse, solution by iterative methods has been studied extensively [1, 5–16]. Much attention has focused on the Navier-Stokes problem; see, for example, [1, 2, 5, 17, 18]. The techniques for solving systems like (1) are so numerous that it is almost impossible to give an overview. In addition to the methods developed specifically for Navier-Stokes problems, existing techniques also include splitting schemes [2, 6, 12, 19, 20], constraint preconditioning [2, 20, 21], Uzawa-type algorithms [2, 12, 22], and (preconditioned) Krylov subspace methods based on (approximations to) the Schur complement [2, 16, 23].

We start with augmentation block triangular preconditioners for the general system (1); see Section 2 for our assumptions. When 𝐴 is nonsingular, results for the general system have been obtained before; for example, Murphy et al. [23, 24] propose the block diagonal Schur complement preconditioner and the block triangular Schur complement preconditioner as follows: ℳ=𝐴00𝐶𝐴⁻¹𝐵^𝑇,ℳ=𝐴0𝐶𝐶𝐴⁻¹𝐵^𝑇. (2) If defined, it has been shown that the preconditioned matrices (cf. [24]) ℳ⁻¹𝒜=𝐼𝐴⁻¹𝐵^𝑇𝐶𝐴⁻¹𝐵^𝑇−1,𝐶0ℳ⁻¹𝒜=𝐼𝐴⁻¹𝐵^𝑇0−𝐼 (3) are diagonalizable and have only three distinct eigenvalues 1, (1±√5)/2 and two distinct eigenvalues 1, −1, respectively.

However, when 𝐴 is singular, it cannot be inverted and the Schur complement does not exist. For symmetric saddle point systems, that is, 𝐶=𝐵 and 𝐴 is symmetric, one possible way of dealing with the systems is by augmentation, that is, by replacing 𝐴 with 𝐴+𝐵^𝑇𝑊⁻¹𝐵, where 𝑊 is an 𝑚×𝑚 symmetric positive definite weight matrix [4, 7, 9, 14, 17, 18, 25–27]. Recently, for symmetric saddle point systems with (1,1) block that has a high nullity, Greif and Schötzau [25, 26] studied the application of the following block diagonal preconditioner used with the MINRES solver for the nonsymmetric saddle point systems (1): ℳ=𝐴+𝐵^𝑇𝑊⁻¹.𝐵00𝑊 (4) They have shown that if the nullity of 𝐴 is 𝑚, which is the highest possible nullity, then the preconditioned matrix ℳ⁻¹𝒜 has only two distinct eigenvalues 1 and −1. Thus, a preconditioned minimal residual Krylov subspace iterative method such as MINRES converges within two iterations.

Recently, Rees and Greif [4] presented a block triangular preconditioner as follows: 𝒯=𝐴+𝐵^𝑇𝑊⁻¹𝐵𝑘𝐵^𝑇0𝑊 (5) and has shown that if 𝑘=−1, the preconditioned matrices 𝒯⁻¹𝒜 have 𝜆=1 with multiplicity 𝑛−𝑚 and (1±√5)/2 each with multiplicity 𝑝. The remaining eigenvalues lie in the intervals [(1−√5)/2,0]⋃[1,(1+√5)/2]. For general nonsymmetric saddle point problems with nonsingular (1,1) block, some block structured preconditioning approaches are also available [3, 4, 18, 22, 28–33].

In this paper, we propose two new block triangular preconditioners 𝓐=𝐴+𝐵^𝑇𝑊⁻¹𝐶𝑘𝐵^𝑇,0−𝑊𝓐_𝑘=𝐴+𝐵^𝑇𝑊⁻¹𝐶𝑘𝐵^𝑇0(1−𝑘)𝑊,𝑘≠1, (6)(7) where 𝑘 is a scalar, and 𝑊 is an 𝑚×𝑚 symmetric positive definite weight matrix.

The remainder of this paper is organized as follows. In Section 2, we discuss two block triangular preconditioners for solving nonsymmetric saddle point systems, and algebraic properties are studied too. In Section 3, numerical experiments are provided to validate our analysis in Section 2. In Section 4, we draw some conclusions.

2. Block Triangular Preconditioners

We will adopt the general notation 𝒜=𝐴𝐵^𝑇𝐶0 (8) to represent the nonsymmetric saddle point matrix of (1). We assume that 𝐴 is symmetric and positive semidefinite with nullity 𝑠 and that 𝐵 is of size 𝑚×𝑛 and has full row rank. Note that the assumption that 𝒜 is nonsingular implies that rank(𝐵)=rank(𝐶)=𝑚, which we use in our analysis below.

We next give two lemmas in [2, 27].

Lemma 1 (see [2, 27]). The nonsymmetric saddle point matrix (7) 𝒜=𝐴𝐵^𝑇𝐶0 (9) is nonsingular if and only if the following conditions are satisfied: rank(𝐵)=rank(𝐶)=𝑚,𝑉^𝑇_𝐵2𝐴𝑉_𝐶2𝑛𝑜𝑛𝑠𝑖𝑛𝑔𝑢𝑙𝑎𝑟, (10) where 𝑉_𝐵2 and 𝑉_𝐶2 are bases of 𝒩(𝐵) and 𝒩(𝐶), respectively. Obviously, the nonsingularity of 𝑉^𝑇_𝐵2𝐴𝑉_𝐶2 implies 𝒩(𝐴)(𝐶)𝒩𝐴∩𝒩=0,_𝑇∩𝒩(𝐵)=0. (11)

Lemma 2 (see [2, 27]). If the saddle point-type matrix 𝒜 in (7) is nonsingular, then the rank of the matrix 𝐴 is at least 𝑛−𝑚, and hence its nullity is at most 𝑚.

2.1. Block Triangular Preconditioners 𝓐

We first consider the preconditioner 𝓐=𝐴+𝐵^𝑇𝑊⁻¹𝐶𝑘𝐵^𝑇.0−𝑊 (12) It is easy to see that the eigenvalues of the preconditioned matrix 𝓐⁻¹𝒜 satisfy the generalized eigenvalue problem 𝐴𝐵^𝑇𝑢𝐶0𝑣=𝜆𝐴+𝐵^𝑇𝑊⁻¹𝐶𝑘𝐵^𝑇𝑢𝑣.0−𝑊 (13) The second block row gives 𝑣=−(1/𝜆)𝑊⁻¹𝐶𝑢, and substituting it into the first block equation gives 𝜆(𝜆−1)𝐴𝑢+𝜆²−𝜆𝑘+1𝐵^𝑇𝑊⁻¹𝐶𝑢=0. (14) Regardless of the choice of 𝑘, we see that 𝜆=1 with algebraic multiplicity 𝑚−𝑚. From the nullity of 𝐴 it follows that there are 𝑠 linearly independent null vectors of 𝐴. For each such null vector, we can find two 𝜆 values satisfying 𝜆²−𝜆𝑘+1=0. Thus, we have 𝜆_±=𝑘±𝑘²−42, (15) each with algebraic multiplicity 𝑠. The remaining 2(𝑚−𝑠) eigenvalues satisfy as following: −𝜆²−𝜆𝜆²−𝜆𝑘+1𝐴𝑢=𝐵^𝑇𝑊⁻¹𝐶𝑢. (16) Therefore, we rewrite (16) as 𝐵^𝑇𝑊⁻¹𝐶𝑢=𝜇𝐴𝑥, (17) where −(𝜆²−𝜆)/(𝜆²−𝜆𝑘+1)=𝜇. Thus, we have ±𝜆=1+𝑘𝜇1+𝑘𝜇²−4𝜇1+𝜇2.1+𝜇 (18) Since 𝐵^𝑇𝑊⁻¹𝐶+𝐴 is assumed to be nonsingular, the matrix pencil 𝐵^𝑇𝑊⁻¹𝐶+𝜀𝐴 is regular (cf. [27]). This expression gives an explicit formula in terms of the generalized eigenvalues of (16) and can be used to identify the intervals in which the eigenvalues lie.

To illustrate this, we consider the case 𝑘=1, which corresponds to setting the (1,2) block of the preconditioner to be 𝐴^𝑇. The preconditioned matrix 𝓐⁻¹𝒜 has 𝜆=1 with multiplicity 𝑛−𝑚, and (1±√3𝑖)/2 each with multiplicity 𝑝. By (18) we have 1𝜆=2±−3𝜇²−2𝜇+12.1+𝜇 (19) It is worth noting that since 𝐴 is typically highly singular, many of the generalized eigenvalues are large, in which case the corresponding eigenvalues 𝜆_±=(1±√3𝑖)/2 are bounded away from zero. Since the 2𝑠 eigenvalues 𝜆_±=(𝑘±𝑘²−4)/2 are unbounded as 𝑘 goes to ∝, we conclude that 𝑘 should be of moderate size.

2.2. Block Triangular Preconditioner 𝓐₁

We next consider the preconditioner 𝓐_𝑘=𝐴+𝐵^𝑇𝑊⁻¹𝐶𝑘𝐵^𝑇0(1−𝑘)𝑊,𝑘≠1, (20) where 𝑘≠1 is a scale, 𝑊 is an 𝑚×𝑚 symmetric positive definite weight matrix, and 𝐴+𝐵^𝑇𝑊⁻¹𝐶 is nonsingular.

The following theorem provides the spectrum results of the preconditioned matrix 𝓐_𝑘⁻¹𝒜.

Theorem 3. Let 𝒜 be nonsingular and let its (1,1) block 𝐴 be singular with nullity 𝑠. Then 𝜆=1 is an eigenvalue of 𝓐_𝑘⁻¹𝒜 of geometric multiplicity 𝑛, and 𝜆=1/(𝑘−1) is an eigenvalue of geometric multiplicity 𝑠. The remaining 𝑚−𝑠 eigenvalues satisfy the relation 𝜇𝜆=−(,1−𝑘)𝜇+1 (21) where 𝜇 are some 𝑚−𝑠 generalized eigenvalues of the following generalized eigenvalue problem: 𝐵^𝑇𝑊⁻¹𝐶𝑥=𝜇𝐴𝑥. (22)

Proof. Suppose that 𝜆 is an eigenvalue of 𝓐_𝑘⁻¹𝒜, whose eigenvector is (𝑢^𝑇,𝑣^𝑇)^𝑇. Then we have 𝓐_𝑘⁻¹𝒜𝑢𝑣𝑢𝑣=𝜆. (23) Furthermore, it satisfies the generalized eigenvalue problem 𝐴𝐵^𝑇𝑢𝐶0𝑣=𝜆𝐴+𝐵^𝑇𝑊⁻¹𝐶𝑘𝐵^𝑇0()𝑊𝑢1−𝑘𝑣 (24) or 𝐴𝑢+𝐵^𝑇𝑣=𝜆𝐴+𝐵^𝑇𝑊⁻¹𝐶𝑢+𝜆𝑘𝐵^𝑇()𝑣,𝐶𝑢=𝜆1−𝑘𝑊𝑣. (25) As 𝓐_𝑘 is nonsingular, 𝜆≠0. The second equality gives that 𝑣=(1/𝜆(1−𝑘))𝑊⁻¹𝐶𝑢, and substituting it into the first equality gives (𝜆−1)𝜆(1−𝑘)𝐴𝑢−(𝜆𝑘−𝜆−1)𝐵^𝑇𝑊⁻¹𝐶𝑢=0. (26) It is straightforward to see that any vector 𝑢∈𝑅^𝑛 satisfies (26) with 𝜆=1, and thus 𝜆=1 is an eigenvalue of 𝓐_𝑘⁻¹𝒜. We claim that the eigenvalue 𝜆=1 has geometric multiplicity 𝑛.
If 𝜆≠1, then from (26) we have 𝜆(1−𝑘)𝐴𝑢+(𝜆𝑘−𝜆−1)𝐵^𝑇𝑊⁻¹𝐶𝑢=0. (27) If 𝑢∈𝒩(𝐴), then from (27) we have (𝜆𝑘−𝜆−1)𝐵^𝑇𝑊⁻¹𝐶𝑢=0. (28) From Lemma 1 and (28), 𝐴 is singular with nullity 𝑠, we know that 𝜆=1/(𝑘−1) is an eigenvalue of 𝓐_𝑘⁻¹𝒜 of geometric multiplicity 𝑠.
We have determined 𝑛+𝑠 eigenvalues. Now we consider the remaining 𝑚−𝑠 eigenvalues of 𝓐_𝑘⁻¹𝒜.
Suppose that 𝜆≠1,1/(𝑘−1), then from (27) we have 𝐵^𝑇𝑊⁻¹𝐶𝑢=𝜇𝐴𝑢, (29) where 𝜇=𝜆(1−𝑘)/(𝜆𝑘−𝜆−1), which implies that 𝜆=−𝜇/(1−𝑘)(𝜇+1). Since 𝐵^𝑇𝑊⁻¹𝐶+𝐴 is assumed to be nonsingular, the matrix pencil 𝐵^𝑇𝑊⁻¹𝐶+𝜀𝐴 is regular (cf. [27]). Thus, the generalized eigenvalue problem (29) is well posed.

If 𝑘=0, then we have the following corollary.

Corollary 4. Let 𝒜 be nonsingular and let its (1,1) block 𝐴 be singular with nullity 𝑠. Then 𝜆=1 is an eigenvalue of 𝓐₀⁻¹𝒜 of geometric multiplicity 𝑛, and 𝜆=−1 is an eigenvalue of geometric multiplicity 𝑠. The remaining 𝑚−𝑠 eigenvalues satisfy the relation 𝜇𝜆=−,𝜇+1 (30) where 𝜇 are some 𝑚−𝑠 generalized eigenvalues of the following generalized eigenvalue problem: 𝐵^𝑇𝑊⁻¹𝐶𝑥=𝜇𝐴𝑥. (31)

From Theorem 3 and Corollary 4 we know that the higher the nullity of 𝐴 is, the stronger the clustered eigenvalues of 𝓐₀⁻¹𝒜 are. When the nullity of 𝐴 is 𝑚, its at most value (cf. Lemma 2), we have the following result.

Corollary 5. Let 𝒜 be nonsingular and let its (1,1) block 𝐴 be singular with nullity 𝑚. Then the preconditioned matrix 𝓐₀⁻¹𝒜 is diagonalizable and has precisely two eigenvalues 𝜆=1 of geometric multiplicity 𝑛 and 𝜆=−1 of geometric multiplicity 𝑚.

If 𝑘=2, we have the following corollary.

Corollary 6. Let 𝒜 be nonsingular and let its (1,1) block 𝐴 be singular with nullity 𝑠. Then 𝜆=1 is an eigenvalue of 𝓐₂⁻¹𝒜 of geometric multiplicity 𝑛+𝑠. The remaining 𝑚−𝑠 eigenvalues satisfy the relation 𝜇𝜆=,𝜇+1 (32) where 𝜇 are some 𝑚−𝑠 generalized eigenvalues of the following generalized eigenvalue problem: 𝐵^𝑇𝑊⁻¹𝐶𝑥=𝜇𝐴𝑥. (33)

From Theorem 3 and Corollary 6, if the nullity of 𝐴 is 𝑚, we have the following corollary.

Corollary 7. Let 𝒜 be nonsingular and let its (1,1) block 𝐴 be singular with nullity 𝑚. Then the preconditioned matrix 𝓐₂⁻¹𝒜 has only precisely one eigenvalue 𝜆=1 of geometric multiplicity 𝑛+𝑚.

Remark 8. From Corollaries 5 and 7, if the nullity of 𝐴 is 𝑚, then it is readily seen that the preconditioned matrix 𝓐₀⁻¹𝒜 has precisely two distinct eigenvalues and that the preconditioned matrix 𝓐₂⁻¹𝒜 has only precisely one eigenvalue. Thus, we know that any preconditioned Krylov subspace method such as GMRES terminates in at most two steps if roundoff errors are ignored.

3. Numerical Experiments

In this section we present numerical experiments to illustrate the performance of the two preconditioners when they are implemented exactly or approximately.

For our experiments, as in [19, 27] we construct the saddle point-type matrix 𝒜 from reforming a matrix 𝒜 of the following form: 𝒜=𝐴𝐵^𝑇,𝐵0 (34) where the matrix 𝒜 arises from the discretization of the Stokes problem as follows: −𝜇Δ𝐮+∇𝑤=𝑓,inΩ,∇⋅𝐮=𝑔,inΩ,𝐮=0,on𝜕Ω,_Ω𝑤(𝑥)𝑑𝑥=0, (35) where Ω=(0,1)×(0,1)⊂𝑅², 𝜕Ω is the boundary of Ω, Δ is the componentwise Laplace operator, 𝑢 is a vector-valued function representing the velocity, and 𝑤 is a scalar function representing the pressure (cf [19]). Then the matrices 𝐵 and 𝐵^𝑇 are replaced by a random matrix 𝐶 with the same sparsity as 𝐵 and a random matrix 𝐵^𝑇 with the same sparsity as 𝐵^𝑇, respectively. Furthermore, 𝐶(1∶𝑚,1∶𝑚) and 𝐵(1∶𝑚,1∶𝑚) are, respectively, replaced by 𝐶₁=𝐶(1∶𝑚,1∶𝑚)−(3/2)𝐼_𝑚 and 𝐵₁=𝐵(1∶𝑚,1∶𝑚)−(3/2)𝐼_𝑚, such that 𝐶₁ and 𝐵₁ are nonsingular. Denote 𝐶₂=𝐶(1∶𝑚,𝑚+1∶𝑛) and 𝐵₂=𝐵(1∶𝑚,𝑚+1∶𝑛); then we have 𝐶=[𝐶₁,𝐶₂] and 𝐵=[𝐵₁,𝐵₂] with 𝐵₁,𝐶₁∈𝑅^𝑚,𝑚 and 𝐵₂,𝐶₂∈𝑅^{𝑚,𝑛−𝑚}. Obviously, the resulting saddle point-type matrix 𝒜=𝐴𝐵^𝑇𝐶0 (36) satisfies rank(𝐶)=rank(𝐵)=𝑚. (37) From the matrix 𝒜 in (36) we construct the following four saddle point-type matrices: 𝒜_𝑖=𝐴_𝑖𝐵^𝑇𝐶0,𝑖=1,…,4, (38) where 𝐴_𝑖 is constructed from 𝐴 by making its 𝑖×(𝑚/4) rows and columns with zero entries. Note that 𝐴_𝑖 is semipositive real and its nullity is 𝑖×(𝑚/4).

In the first numerical example, the matrix 𝑊 in the augmentation block preconditioners is taken as 𝑊=(1/𝛾)𝐼_𝑚, whereas the positive parameter 𝛾 is taken as (cf. [17]) 𝐴𝛾=^𝑖1‖𝐵‖₁‖𝐶‖₁. (39)

Theorem 3 and Corollaries 4–7 guarantee increased spectral clustering of the preconditioned matrix 𝓐_𝑘⁻¹𝒜 when 𝑘=2. As expected by the previous theorems and corollaries, Figures 1 and 2 plot, respectively, the two preconditioned matrices, each figure includes four subfigures that plot, respectively, the eigenvalues of four preconditioned matrices resulting from four matrices 𝒜_𝑖, 𝑖=1,…,4, operated upon by exactly the same preconditioner.

(a)

(b)

(c)

(d)

Figure 1 

The eigenvalue distribution of preconditioned systems 𝒯−1𝒜 with 𝑘=2 and 𝑚+𝑛=2700.

(a)

(b)

(c)

(d)

Figure 2 

The eigenvalue distribution of preconditioned systems 𝓐−1𝒜 with 𝑘=2 and 𝑚+𝑛=2700.

From these figures we can clearly see that the higher the nullity of the (1,1) block, the more strongly the eigenvalues of the preconditioned matrices are clustered. Furthermore, we note that Figures 1 and 2 show that 𝜆=1 has very high multiplicity when the nullity of the (1,1) block 𝐴_𝑖 is 𝑚.

Next, we study the iteration numbers and iteration time for 𝒜₃ and 𝒜₄. Table 1 shows results for applying GMRES of block triangular preconditioners 𝓐. From Table 1, we can see that the preconditioned GMRES of 𝓐⁻¹𝒜 is more efficient when 𝑘=2, and iteration numbers are slightly changed by the change of parameter 𝑘 when the nullity of 𝐴 is 𝑚 (see Tables 1 and 2).

Table 1 

Iteration time and iteration numbers of GMRES of block triangular preconditioner 𝓐 for different 𝑘 with 𝑚+𝑛=1200.

Table 2 

Iteration time and iteration numbers of GMRES of preconditioned matrices 𝒯⁻¹𝒜 and 𝓐⁻¹𝒜 for different 𝑘=2 with nullity =3𝑚/4.

We next will use preconditioned GMRES to solve, respectively, the following four saddle point-type systems: 𝒜_𝑖𝑥=𝑓_𝑖,𝑖=3,4, (40) where the right-hand side 𝑓_𝑖 is taken such that the exact solution 𝑥 is equal to (1,1…,1)^𝑇. The stopping criterion is ‖𝑟^(𝑘)‖/‖𝑓_𝑖‖≤10⁻⁶ and 𝑟^(𝑘) is the residual vector after 𝑘th iteration. In this section, let 𝑘=2; we give some comparison results of the iteration numbers and iteration time for five preconditioners (see Tables 2 and 3). From Table 2, when the nullity of 𝐴 is 3𝑚/4, we see that the preconditioned GMRES with the preconditioner 𝓐 is more efficient than that of the preconditioner 𝒯, but from Table 3, when the nullity of 𝐴 is 𝑚, we see that the preconditioned GMRES with the preconditioner 𝓐 is essentially similar to the preconditioner 𝒯 in iteration time and iteration numbers.

Table 3 

Iteration time and iteration numbers of GMRES of preconditioned matrices 𝒯⁻¹𝒜 and 𝓐⁻¹𝒜 for different 𝑘=2 with nullity =𝑚.

All the numerical experiments were performed with MATLAB 7.0. The machine we have used is a PC-AMD, CPU T7400 2.2 GHz process.

4. Conclusions

In this paper, we have analyzed the spectral properties as well as the computational performance of two types of block triangular augmentation preconditioners for saddle point-type matrices with a highly singular (1,1) block. Complete theoretical analysis shows that all eigenvalues of the preconditioned matrices are strongly clustered. A good parameter choice may substantially reduce the iteration numbers and iteration time. In particulare, we have shown that in cases where the (1,1) block has high nullity, convergence for each of the two preconditioned GMRES iterative methods is guaranteed to be almost immediate. Numerical experiments are also reported for illustrating the efficiency of the presented preconditioners. From the views of theories and applications, the presented preconditioners are better than the previous results.

Acknowledgments

The author wishes to thank the editor Professor Kamal Mostafa and the anonymous referees for their helpful suggestions which improved this paper. This research was supported by the National Natural Science Foundation of China (no. 11071079), the Ningbo Natural Science Foundation (2012A610037) and the Fundamental Research Funds for the Central Universities.