Spectral Analysis, Properties and Nonsingular Preconditioners for Singular Saddle Point Problems

1 Introduction

In this work we shall mainly study the spectral behavior and preconditioners of the generalized singular or nonsingular saddle point problems of the form

where A∈ℝn×n is a symmetric positive definite matrix, C∈ℝm×m is a symmetric positive semi-definite matrix, and B is a general matrix in ℝm×n, with m≤n (possibly m≪n); f∈ℝn and g∈ℝm are two given vectors.

Systems of the form (1.1) arise from a variety of scientific and engineering applications, such as computational fluid dynamics, optimal control problems, constrained optimizations, weighted least-squares formulations, mixed or hybrid finite element approximations of partial differential equations, electronic networks and computer graphics, boundary element discretization based on domain decomposition, and so on; see, e.g., [14, 11, 29, 20, 23, 21, 10, 22, 6] and the references therein.

Although the linear system (1.1) can be equivalently and easily transformed into a symmetric linear system, the coefficient matrix of non-symmetric form in (1.1) is semi-positive real and positive semi-stable (see [5] for more details). So the use of the non-symmetric form (1.1) can be advantageous when using certain Krylov subspace methods, like restarted GMRES [5], and has motivated many studies on the numerical solution of the generalized saddle point systems of the form (1.1), see, e.g., [18, 7, 8, 2, 5, 3, 31, 32, 36]. This will be also the central focus of the current work, namely, we shall mainly study the spectral properties of the generalized non-symmetric saddle point system (1.1), especially for the case that it is singular, and propose some effective nonsingular preconditioners for the singular systems and analyze the spectral behavior of the preconditioned systems. When the generalized saddle point systems of the form (1.1) are large and sparse, iterative methods appear to be more attractive than direct methods. However, most existing iterative methods, including Krylov subspace-type methods, converge quite slowly, or even fail to converge if they are ineffectively preconditioned. Therefore, great attention has been paid to the preconditioned iterative methods and the analysis on their convergence behavior in the last decade; see, e.g., [16, 4, 30, 17, 25, 26, 27]. Quite many preconditioners have been proposed and studied recently for the saddle point systems of the form (1.1), such as symmetric indefinite preconditioners [7, 8], inexact constraint preconditioners [8, 28, 3, 36], block triangular preconditioners [13, 34, 21, 33], shift-splitting preconditioners [5, 15, 35] and primal-based penalty preconditioners [19]. In particular, when the saddle point system (1.1) is generated by some special discretizations, e.g., finite element or boundary element methods based on non-overlapping domain decomposition, effective preconditioners can be constructed by making full use of the special structure of discretizations [14, 9]. But we emphasize that most of the above studies were done for the cases when the saddle point systems (1.1) are nonsingular, and few analytical results on preconditioning effects are available for the cases when the saddle point systems (1.1) are singular, which have wide applications in many areas as we mentioned earlier. This will be the central task of our current work, and several new analysis tools and techniques are introduced for the analytical purpose. As it is widely known, an effective and efficient preconditioner is mostly expected to have the clustered eigenvalues for its resultant preconditioned system. In order to see such clustering behavior, one needs some careful and sharp estimates of the spectral bounds for the preconditioned system. For non-symmetric generalized saddle point matrices of the form

in (1.1), there are already quite many investigations available in literature; see [7, 8, 32, 34, 2, 3, 24]. Interesting estimates were established in [7] on spectral bounds of the complex eigenvalues of 𝒜, and further improved in [32] from the previous lower bound (ρ1+μ1)/2 to the lower bound max⁡{(ρ1+μ1)/2,μ1-δm}, where ρ1 and μ1 are the minimum eigenvalues of A and C, respectively, and δm is the maximum eigenvalue of B⁢BT. But the matrix C is often singular or even vanishes in many applications, so we have μ1=0 and the improved estimates of [32] reduce to the original ones in [7]. The first results of this work are to improve the estimates in [32] for the important cases when matrix C is singular or zero. With the help of those improved results, we shall derive some explicit and sharp bounds on the spectra of the generalized non-symmetric singular saddle point matrices of the form (1.1). Then we propose new nonsingular preconditioners for solving the generalized singular saddle point problems, and establish the spectral bounds of the preconditioned systems.

The arrangement of this work is in the following order. We first present some improved estimates of the spectra for the generalized singular or nonsingular saddle point matrix 𝒜 in Section 2, then study more specific spectral properties of singular matrices 𝒜 in Section 3. We will then propose a new nonsingular preconditioner for the singular saddle point system in (1.1) in Section 4, and derive the spectral bounds of the preconditioned systems in Section 5. Finally, in Section 6, the nonsingular preconditioners are applied to solve the singular finite element saddle point systems arising from the discretization of the Stokes problems to demonstrate their effectiveness and efficiencies.

2 Estimates of Complex Eigenvalues of Singular or Nonsingular Saddle Point Matrices

In this section, we make an effort to improve the following estimates of eigenvalues of matrix 𝒜 in [32], which will be important to our subsequent studies of the spectral distributions of matrix 𝒜 and its preconditioned system.

It is often encountered in many applications that the matrix C is singular. In this case, Theorem 2.1 implies that the lower bound of the real part of any complex eigenvalue of 𝒜 reduces to (ρ1+μ1)/2, which deteriorates to the same result in [7]. In the next theorem we present some estimates that improve the ones in Theorem 2.1 for the important case that the matrix C is singular. There are no existing analytical tools to apply for this singular case. Instead, we choose to consider the eigensystem of the saddle point equation (1.1) directly, from where we solve the variable x in terms of the variable y and substitute into the second equation of the eigensystem. Then we separate the real and imaginary part of the quadratic form of the resultant equation, and achieve the desired estimates by using the important relation between the real and imaginary parts.

Since A is symmetric and positive definite, we have ρ1>0. Theorem 2.2 then shows that our new estimates improve the previous estimates in Theorem 2.1 for the cases that C is singular, provided ρ1-δm>0. The combination of Theorem 2.2 and Theorem 2.1 leads to the following estimates.

3 Spectral Behavior of Singular Saddle Point Matrices

In recent years, several studies have been carried out to analyze the spectral properties of the saddle point matrix 𝒜 of form (1.2) when 𝒜 is nonsingular; see, e.g., [7, 8, 1]. But it is often encountered in many applications that matrix 𝒜 is singular. To the best of our knowledge, not much effort has been made yet in literature to study this important singular case. This is our major focus of this section. We shall first present some estimates of the rank of the singular saddle point matrix 𝒜, which help to determine the multiplicity of the eigenvalue 0. Then we will establish some important estimates of both real and complex parts of all eigenvalues of the singular matrix 𝒜. Those estimates are the essential tools for our analysis in Section 4 on the spectral behavior of the preconditioned systems for singular saddle point matrices by nonsingular preconditioners.

By Theorem 3.1, we know that rank⁡[𝒜]=n+m if rank⁡[C]=m or rank⁡[B]=m. This implies immediately that rank⁡[C]<m, μ1=0 and rank⁡[B]<m if 𝒜 is singular. This leads to the following estimate of the multiplicity of the eigenvalue 0 of 𝒜.

For the real eigenvalues of the singular saddle point matrix 𝒜, we can get the following results by a similar argument as that used in [1, Proposition 2.2] and [31, Theorem 2.1].

The following results are direct consequences of Theorems 3.2 and 3.3.

If C=0, then Corollary 3.4 can be much simplified as stated in the following corollary.

4 Nonsingular Preconditioners for Singular Saddle Point Matrices

Intensive efforts have been made for the studies of preconditioners for nonsingular saddle point systems of form (1.1); see, e.g., [7, 8, 28, 3, 13, 34, 21, 33, 5, 35, 19]. On the other hand, few results are available for the analysis on preconditioners for singular saddle point systems, and we are only aware of the results in [36], which considered the simple case with C=0 in (1.1). Moreover, only singular preconditioners were analyzed in [36] for the singular saddle point system (1.1), but which are not always reliable and reasonable, as the preconditioned systems may not be equivalent to the original system (1.1) if singular preconditioners are used. In this section, we present and study some nonsingular preconditioners 𝒦 for the saddle point matrix 𝒜 in (1.1) when 𝒜 is singular. We shall consider the preconditioners of the general form

where PA and PC are symmetric positive definite approximations of A and C, respectively. Then we shall solve the following preconditioned system:

or

where z=(xT,yT)T and b=(fT,-gT)T. Since

we need only to solve two subsystems PA⁢x=rx and (PC+B⁢PA-1⁢BT)⁢y=ry when applying the GMRES method coupled with the preconditioner 𝒦 to the saddle point problem (1.1). In addition, we may easily see that B is not of full row rank when 𝒜 is singular, so it is natural to require a symmetric positive definite PC in (4.1). By Theorem 3.1, we know that 𝒦 is nonsingular. Some simple choices of PA and PC can be, e.g., PA=diag⁡(A) and PC=α⁢I with a positive parameter α.

The preconditioners of form (4.1) were considered in [8] for the case when 𝒜 is nonsingular, and PC was taken to be the special matrix PC=B⁢PA-1⁢BT-PS, where PS is a symmetric positive definite approximation of the Schur complement S. But we would emphasize that the singular case considered here is much more technical than the nonsingular system. One main technical difficulty one may encounter in solving a singular system by GMRES is that the iteration process may break down at some step. Some conditions were derived in [12] under which the GMRES iterates converge safely to a solution of a singular system. For our subsequent study of the convergence of the GMRES method for solving the preconditioned system (4.2), we introduce the following convergence result from [12, Theorem 2.4].

Using the above convergence theory, we can derive the following result when GMRES is applied for solving the singular saddle point problem (1.1).

Similarly to the results in Theorem 4.2, we have the following convergence for GMRES applied for the preconditioned system (4.3).

Although Theorems 4.2 and 4.4 guarantee the convergence of GMRES when applied to the preconditioned system (4.2) and (4.3), respectively, its convergence rate depends strongly on the spectral properties of the preconditioned matrix. It will be our central task in the subsequent section to study the spectral properties of the preconditioned matrix. For these studies, we remark that the general case C≠0 considered here is much more technical than the simple case C=0 studied in [36], and the analysis of [36] clearly does not apply to the current case with C≠0. Moreover, unlike the results in [36], which gave only the estimates of the real part of the eigenvalues of the preconditioned system in terms of the spectrum of the preconditioned system for A in (1.1), we shall derive the lower and upper bounds of both the real and complex eigenvalues of the preconditioned matrix.

5 Spectral Estimates of the Preconditioned Matrix

There are two major analysis techniques that enable the study of the spectral behavior of the preconditioned matrix 𝒜⁢𝒦-1 in (4.3). First, we manage to work out two effective similar transformations, which transform our preconditioned system 𝒜⁢𝒦-1 into a matrix, see ℳ in (5.2) below, that has exactly the same block form as the original saddle point matrix 𝒜 in (1.2). Then we can apply our results of Section 3 for the saddle point matrix 𝒜 to derive the lower and upper bounds of the real and complex eigenvalues of the preconditioned system 𝒜⁢𝒦-1. But this application is not straightforward. We thus have to establish the spectral estimates for each of the matrices Q⁢QT, H and the Schur complement W=:H+QAp-1QT (see (5.2)–(5.4) below) associated with the new saddle point matrix ℳ in (5.2), and several delicate new techniques are introduced here.

In the remainder of this section, we will investigate the spectral properties of 𝒜⁢𝒦-1 by using the results we developed in Section 3. To this end, we first define a matrix

We can directly verify that 𝒜⁢𝒦-1 is similar to (𝒫-1⁢𝒜⁢𝒫-1)⁢(𝒫⁢𝒦-1⁢𝒫), i.e.,