Inexact modified positive-definite and skew-Hermitian splitting preconditioners for generalized saddle point problems

Abstract

In this article, to better implement the modified positive-definite and skew-Hermitian splitting preconditioners studied recently (Numer. Algor., 72 (2016) 243–258) for generalized saddle point problems, a class of inexact modified positive-definite and skew-Hermitian splitting preconditioners is proposed with improved computing efficiency. Some spectral properties, including the eigenvalue distribution, the eigenvector distribution, and an upper bound of the degree of the minimal polynomial of the inexact modified positive-definite and skew-Hermitian splitting preconditioned matrices are studied. In addition, a theoretical optimal inexact modified positive-definite and skew-Hermitian splitting preconditioner is obtained. Numerical experiments arising from a model steady incompressible Navier–Stokes problem are used to validate the theoretical results and illustrate the effectiveness of this new class of proposed preconditioners.

Introduction

In the solution of large, sparse generalized saddle point problems, finding efficient preconditioners is crucial to obtaining an iterative solution. Consider a matrix of the form

$$\mathcal{A}\mathrm{x}\equiv \left[\begin{array}{cc}A\hfill & B^{*}\hfill \\ \multicolumn{1}{c}{-B}& C\hfill \end{array}\right]\left[\begin{array}{c}x\hfill \\ \multicolumn{1}{c}{y}\end{array}\right]=\left[\begin{array}{c}f\hfill \\ \multicolumn{1}{c}{g}\end{array}\right]\equiv b$$

(1)

where $A\in {\mathbb{C}}^{n\times n}$ is a non-Hermitian positive-definite matrix (i.e. its Hermitian part $A+A^{*}$ is positive definite), $B\in {\mathbb{C}}^{m\times n}$ is a rectangular matrix with full row rank, $C\in {\mathbb{C}}^{m\times m}$ is Hermitian positive-semidefinite and $m\le n$ . Here, $(·{)}^{*}$ denotes the conjugate transpose of the corresponding matrix or vector. From the algebraic properties studied in Benzi et al., ¹ we know that the non-Hermitian positive-semidefinite generalized saddle point matrix $\mathcal{A}$ is nonsingular. Thus, the linear system (1) has a unique solution.

Linear systems of the form (1) arise in a number of engineering applications, including piezoelectric structure computation, mixed finite element solutions of the Navier–Stokes equations, and constrained optimization problems (see Benzi et al. ¹ for a general discussion). Since the coefficient matrix $\mathcal{A}$ is large, sparse, nonsymmetric, and often ill-conditioned, many traditional iterative methods converge to a solution very slowly. Thus, preconditioning is in most cases indispensable for obtaining an iterative solution of equation (1). ² In the past few decades, many preconditioning techniques have been proposed and studied for such generalized saddle point problems (1). For an overview of these preconditioning techniques, see Benzi et al. ¹ and Pestana and Wathen ³ and the references therein.

Let $\alpha >0$ be a real positive constant and I be the identity matrix with suitable dimension. Based on the Hermitian and skew-Hermitian splitting (HSS) preconditioner studied in Bai et al., ⁴ Benzi and Golub, ⁵ and Bai and Golub, ⁶ Pan et al.7 first split the generalized saddle point matrix $\mathcal{A}$ into

$$\mathcal{A}=\mathcal{P}+\mathcal{S}=\left[\begin{array}{cc}A\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& C\hfill \end{array}\right]+\left[\begin{array}{cc}0\hfill & B^{*}\hfill \\ \multicolumn{1}{c}{-B}& 0\hfill \end{array}\right]$$

(2)

where $\mathcal{P}$ is a block diagonal positive-semidefinite matrix and $\mathcal{S}$ is a skew-Hermitian matrix and then proposed a matrix splitting preconditioner as follows

$$\begin{array}{l}{\mathcal{P}}_{DPSS}=\frac{1}{2\alpha }(\alpha I+\mathcal{P})(\alpha I+\mathcal{S})\\ =\frac{1}{2\alpha }\left[\begin{array}{cc}\alpha I+A& 0\\ 0& \alpha I+C\end{array}\right]\left[\begin{array}{cc}\alpha I& B^{\ast }\\ -B& \alpha I\end{array}\right]\end{array}$$

(3)

The above matrix is called the deteriorated positive-definite and skew-Hermitian splitting (DPSS) preconditioner. Evidently, if the (1,1) leading block matrix A is Hermitian, then the DPSS preconditioner is the same as the HSS preconditioner. ⁵ The clustering property of the eigenvalues of the DPSS preconditioned matrix has been studied. ⁷ The unconditional convergence property of the corresponding DPSS iteration method can be found in Shen, ⁸ Bai, ⁹ and Cao et al. ¹⁰ Numerical results show that the DPSS preconditioner performs better than the well-known HSS preconditioner.

Much effort has gone into improving the DPSS preconditioner. For example, based on the relaxed techniques studied in previous works,^11–13 Cao et al.^10,14 proposed a relaxed deteriorated positive-definite and skew-Hermitian splitting (RDPSS) preconditioner for equation (1) when $C=0$ and a simplified Hermitian and skew-Hermitian splitting (SHSS) preconditioner when the submatrices A and C are Hermitian, the resulting preconditioners are closer to the generalized saddle point matrix than the DPSS-based preconditioner. Zhou et al. ¹⁵ proposed a relaxed block triangular splitting preconditioner for generalized saddle point matrix by considering the coefficient matrix as a three-by-three block matrix. As another example, see Cao et al., ¹⁶ a relaxed PSS-based preconditioner was proposed when the submatrix A is non-Hermitian positive definite. Also recently, by introducing another positive parameter $\beta$ , Xie and Ma ¹⁷ proposed a modified positive-definite and skew-Hermitian splitting (MPSS) preconditioner as follows

$$\begin{array}{l}{\mathcal{P}}_{MPSS}=\left[\begin{array}{cc}A& 0\\ 0& I\end{array}\right]\left[\begin{array}{cc}I& \left(A^{-1}+\frac{1}{\alpha }I\right)B^{\ast }\\ -B& \beta I+C\end{array}\right]\\ =\left[\begin{array}{cc}A& B^{\ast }+\frac{1}{\alpha }A{B}^{\ast }\\ -B& \beta I+C\end{array}\right]\end{array}$$

(4)

which can be further decomposed into

$$\begin{array}{l}{\mathcal{P}}_{MPSS}=\left[\begin{array}{cc}A& 0\\ -B& I\end{array}\right]\left[\begin{array}{cc}I& 0\\ 0& \beta I+C+B{A}^{-1}{B}^{\ast }+\frac{1}{\alpha }B{B}^{\ast }\end{array}\right]\\ \left[\begin{array}{cc}I& \left(A^{-1}+\frac{1}{\alpha }I\right)B^{\ast }\\ 0& I\end{array}\right]\end{array}$$

(5)

From equation (4), we can see that the MPSS preconditioner is sufficiently close to the coefficient matrix $\mathcal{A}$ as $\alpha \to +\infty$ and $\beta \to 0$ . Although Xie et al. have proved that the MPSS preconditioned matrix ${\mathcal{P}}_{\mathrm{MPSS}}^{-1}\mathcal{A}$ has a clustered eigenvalue distribution, the factorization (4) indicates that the MPSS preconditioner ${\mathcal{P}}_{\mathrm{MPSS}}$ is very difficult to implement since a sub-system of linear equations with the coefficient matrix $\beta I+C+B{A}^{-1}{B}^{*}+(1/\alpha )B{B}^{*}$ must be solved at each iteration, see the next section for detailed implementation considerations. Just the computation of the coefficient matrix $\beta I+C+B{A}^{-1}{B}^{*}+(1/\alpha )B{B}^{*}$ can be prohibitive. A remedy to such situations is to use inexact variants. It should be noted that the inexact preconditioners are very important in practical computation since their iteration costs are much cheaper than exact preconditioners and the necessary iteration steps do not increase significantly. For the generalized saddle point problem (1), where the submatrix A is Hermitian positive definite, some inexact block triangular preconditioners are proposed in Simoncini, ¹⁸ Cao, ¹⁹ and Jiang et al. ²⁰ and improved eigenvalue bounds of the preconditioned matrices are studied in Cao et al. ²¹ For a general block two-by-two non-Hermitian matrix, some inexact block triangular preconditioners are studied in Bai and Ng, ²² Peng and Li, ²³ and Li et al. ²⁴ For standard saddle point problem where the (1,1) block matrix A is Hermitian positive definite and the (2,2) block matrix C is a zero matrix, an inexact constraint preconditioner is proposed in Bergamaschi et al. ²⁵ and spectral analysis of the inexact constraint preconditioned matrix is studied in Sesana and Simoncini. ²⁶

In this article, we propose a class of inexact modified positive-definite and skew-Hermitian splitting (IMPSS) preconditioners for the generalized saddle point problem (1). From a computational viewpoint, the proposed IMPSS preconditioners are iteratively cheaper than the exact MPSS preconditioners studied in Xie and Ma. ¹⁷ On the theoretical side, the IMPSS preconditioned matrix also exhibits very desirable spectral properties, which allow us to obtain a theoretical optimal modified positive-definite and skew-Hermitian splitting (OIMPSS) preconditioner.

The rest of this article is organized as follows. The IMPSS preconditioner and implementation aspects will be presented in section “The inexact MPSS preconditioner.” In section “Spectral analysis of the IMPSS preconditioned matrix,” we study the spectral properties of the IMPSS preconditioned matrix and introduce the OIMPSS preconditioner. Afterward, numerical experiments arising from a model steady incompressible Navier–Stokes problem are provided in section “Numerical experiments” to illustrate the theoretical results and examine the numerical effectiveness of the new proposed preconditioners.

The inexact MPSS preconditioner

Before proposing a new efficient preconditioner, we first consider the implementation aspects of the MPSS preconditioner. Let $r=[r_1^{*},r_2^{*}{]}^{*}$ ( $r_1\in {\mathbb{C}}^n$ , $r_2\in {\mathbb{C}}^m$ ) be the generalized residual vector, at each iteration step of the MPSS preconditioned iteration method (as in the preconditioned generalized minimal residual [GMRES] iteration method), we need to solve the following generalized residual equation

$${\mathcal{P}}_{\mathrm{MPSS}}z=r$$

(6)

to obtain the current residual $z=[z_1^{*},z_2^{*}{]}^{*}$ ( $z_1\in {\mathbb{C}}^n$ , $z_2\in {\mathbb{C}}^m$ ). From the matrix factorization expression (5) of the MPSS preconditioner, the current residual z can be solved according to the following algorithm.

Algorithm 1

Computation of $z={\mathcal{P}}_{\mathrm{MPSS}}^{-1}r$ :

1. Solve $A{w}_1=r_1$ to obtain $w_1$ ;

2. Solve $(\beta I+C+B{A}^{-1}{B}^{*}+(1/\alpha )B{B}^{*})z_2=B{w}_1+r_2$ to obtain $z_2$ ;

3. Compute $w_2=B^{*}{z}_2$ and solve $A{w}_3=w_2$ to obtain $w_3$ ;

4. Compute $z_1=w_1-(1/\alpha )w_2-w_3$ .

From Algorithm 1, we can see that the implementation of the MPSS preconditioner requires the solution of three sub-linear systems. The primary computational cost is realized in the second step, that is, finding the inverse of the matrix $\beta I+C+B{A}^{-1}{B}^{*}+(1/\alpha )B{B}^{*}$ exactly. To avoid this, an algorithm variant is proposed in Xie and Ma ¹⁷ [Algorithm 3.1], see also the Algorithm 2.

Algorithm 2

Computation of $z={\mathcal{P}}_{\mathrm{MPSS}}^{-1}r$ :

1. Solve $A{w}_1=r_1$ to obtain $w_1$ . Solve $A{q}_i=B^{*}(:,i)$ to obtain $q_i$ ( $i=1,2,\dots ,m$ ). Set $Q=[q_1,q_2,\dots ,q_m]$ , where $B^{*}(:,i)$ denotes the ith column of matrix $B^{*}$ ;

2. Solve $(\beta I+C+\mathrm{BQ}+(1/\alpha )B{B}^{*})z_2=B{w}_1+r_2$ to obtain $z_2$ ;

3. Compute $w_2=B^{*}{z}_2$ and solve $A{w}_3=w_2$ to obtain $w_3$ ;

4. Compute $z_1=w_1-(1/\alpha )w_2-w_3$ .

Comparing Algorithm 2 with Algorithm 1, we can see that an additional sub-linear system with multiple right-hand sides needs to be solved. Besides this, it may also be that the coefficient matrix in the second step may be dense. It is also difficult to solve large sparse problems.

Based on the inexact preconditioners studied in Simoncini, ¹⁸ Bai and Ng, ²² and Bergamaschi et al., ²⁵ by introducing a nonsingular preconditioning matrix $\hat{S}$ , we propose a class of IMPSS preconditioners as follows

$$\begin{array}{l}{\mathcal{P}}_{IMPSS}=\left[\begin{array}{cc}A& 0\\ -B& I\end{array}\right]\left[\begin{array}{cc}I& 0\\ 0& \widehat{S}\end{array}\right]\left[\begin{array}{cc}I& \left(A^{-1}+\frac{1}{\alpha }I\right)B^{*}\\ 0& I\end{array}\right]\\ =\left[\begin{array}{cc}A& B^{*}+\frac{1}{\alpha }A{B}^{*}\\ -B& \widehat{S}-B\left(A^{-1}+\frac{1}{\alpha }I\right)B^{*}\end{array}\right]\end{array}$$

(7)

where we assume that ${\hat{S}}^{-1}$ is easy to obtain. Evidently, if $\hat{S}=\beta I+C+B{A}^{-1}{B}^{*}+(1/\alpha )B{B}^{*}$ , then the IMPSS preconditioner becomes the MPSS preconditioner (5).

The implementation of the IMPSS preconditioner requires the solution of the following generalized residual equation

$${\mathcal{P}}_{\mathrm{IMPSS}}z=r$$

where $z=[z_1^{*},z_2^{*}{]}^{*}$ and $r=[r_1^{*},r_2^{*}{]}^{*}$ ( $z_1,r_1\in {\mathbb{C}}^n$ , $z_2,r_2\in {\mathbb{C}}^m$ ) are the current and the generalized residual vectors, respectively. According to the matrix factorization (7), the above generalized residual equation with respect to the IMPSS preconditioner can be solved according to Algorithm 3.

Algorithm 3

Computation of $z={\mathcal{P}}_{\mathrm{IMPSS}}^{-1}r$ :

1. Solve $A{w}_1=r_1$ to obtain $w_1$ ;

2. Solve $\hat{S}{z}_2=B{w}_1+r_2$ to obtain $z_2$ ;

3. Compute $w_2=B^{*}{z}_2$ and solve $A{w}_3=w_2$ to obtain $w_3$ ;

4. Compute $z_1=w_1-(1/\alpha )w_2-w_3$ .

Comparing Algorithm 3 with Algorithm 1, we can find that only the second step is different. Algorithm 3 may also be easier to implement than Algorithm 2. So, from a computational viewpoint, the proposed IMPSS preconditioner may have much better computing efficiency than the MPSS preconditioner. In section “Numerical experiments,” we will confirm this observation through numerical results.

Spectral analysis of the IMPSS preconditioned matrix

In this section, some spectral properties of the IMPSS preconditioned matrix ${\mathcal{P}}_{\mathrm{IMPSS}}^{-1}\mathcal{A}$ will be studied. Also, a theoretical optimal IMPSS preconditioner is obtained.

Theorem 1

Let $A\in {\mathbb{C}}^{n\times n}$ be a non-Hermitian positive-definite matrix, $B\in {\mathbb{R}}^{m\times n}$ have full row rank, $C\in {\mathbb{R}}^{m\times m}$ be a Hermitian positive-semidefinite matrix. Let $\hat{S}\in {\mathbb{R}}^{m\times m}$ be a nonsingular matrix and $\alpha$ be a positive constant. Let the IMPSS preconditioner ${\mathcal{P}}_{\mathrm{IMPSS}}$ be defined as in equation (7). Then the preconditioned matrix ${\mathcal{P}}_{\mathrm{IMPSS}}^{-1}\mathcal{A}$ has the following:

• An eigenvalue 1 with multiplicity at least n;

• The remaining nonunit eigenvalues are the nonunit eigenvalues of the matrix ${\hat{S}}^{-1}(C+B{A}^{-1}{B}^{*})$ .

In particular, if $\hat{S}=C+B{A}^{-1}{B}^{*}$ , then all eigenvalues of the preconditioned matrix ${\mathcal{P}}_{\mathrm{IMPSS}}^{-1}\mathcal{A}$ are 1.

Proof

From equation (7), by direct computation, we have

$$\begin{array}{c}{\mathcal{P}}_{IMPSS}^{-1}\mathcal{A}=I-{\mathcal{P}}_{IMPSS}^{-1}({\mathcal{P}}_{IMPSS}-\mathcal{A})\\ =I-\left[\begin{array}{cc}I& -\left(A^{-1}+\frac{1}{\alpha }I\right)B^{*}\\ 0& I\end{array}\right]\left[\begin{array}{cc}I& 0\\ 0& {\widehat{S}}^{-1}\end{array}\right]\left[\begin{array}{cc}{A}^{-1}& 0\\ B{A}^{-1}& I\end{array}\right]\left[\begin{array}{cc}0& \frac{1}{\alpha }A{B}^{*}\\ 0& \widehat{S}-C-B\left(A^{-1}+\frac{1}{\alpha }I\right)B^{*}\end{array}\right]\\ =\left[\begin{array}{cccc}I& A^{-1}{B}^{*}-\left(A^{-1}+\frac{1}{\alpha }I\right)B^{*}{\widehat{S}}^{-1}(C+B{A}^{-1}{B}^{*})& 0& {\widehat{S}}^{-1}(C+B{A}^{-1}{B}^{*})\end{array}\right]\end{array}$$

(8)

which clearly shows that the preconditioned matrix ${\mathcal{P}}_{\mathrm{IMPSS}}^{-1}\mathcal{A}$ is a block upper triangular matrix. Thus, the eigenvalue distribution is easily obtained.

From equation (8), we can see that if $\hat{S}=C+B{A}^{-1}{B}^{*}$ , then the IMPSS preconditioned matrix ${\mathcal{P}}_{\mathrm{IMPSS}}^{-1}\mathcal{A}$ becomes the following block upper triangular matrix

$${\mathcal{P}}_{\mathrm{IMPSS}}^{-1}\mathcal{A}=\left[\begin{array}{cc}I\hfill & -\frac{1}{\alpha }{B}^{*}\hfill \\ \multicolumn{1}{c}{0}& I\hfill \end{array}\right]$$

(9)

which immediately shows that all eigenvalues of the preconditioned matrix ${\mathcal{P}}_{\mathrm{IMPSS}}^{-1}\mathcal{A}$ are 1.

From the eigenvalue distribution, we can see that the optimal choice of the preconditioning matrix $\hat{S}$ is the Schur complement matrix $C+B{A}^{-1}{B}^{*}$ . For such a case, we call the resulting matrix the theoretical OIMPSS preconditioner

$$\begin{array}{l}{\mathcal{P}}_{OIMPSS}=\left[\begin{array}{cc}A& 0\\ -B& I\end{array}\right]\left[\begin{array}{cc}I& 0\\ 0& C+B{A}^{-1}{B}^{*}\end{array}\right]\\ \left[\begin{array}{cc}I& \left(A^{-1}+\frac{1}{\alpha }I\right)B^{*}\\ 0& I\end{array}\right]=\left[\begin{array}{cc}A& B^{*}+\frac{1}{\alpha }A{B}^{*}\\ -B& C-\frac{1}{\alpha }B{B}^{*}\end{array}\right]\end{array}$$

(10)

Therefore, to obtain a fast convergence rate of the IMPSS preconditioned Krylov subspace method, we can find the preconditioning matrix $\hat{S}$ that approximates the Schur complement matrix $C+B{A}^{-1}{B}^{*}$ rather than $\beta I+C+B{A}^{-1}{B}^{*}+(1/\alpha )B{B}^{*}$ . From the above theoretical results of the eigenvalue distribution, we can see that the MPSS preconditioner is not the optimal one. Now, we turn to the study of the eigenvector distribution of the IMPSS preconditioned matrix and the OIMPSS preconditioned matrix.

Theorem 2

Let $A\in {\mathbb{C}}^{n\times n}$ be a non-Hermitian positive-definite matrix, $B\in {\mathbb{R}}^{m\times n}$ have full row rank, $C\in {\mathbb{R}}^{m\times m}$ be a Hermitian positive-semidefinite matrix. Let $\hat{S}\in {\mathbb{R}}^{m\times m}$ be a nonsingular matrix and $\alpha$ be a positive constant. Let the IMPSS preconditioner ${\mathcal{P}}_{\mathrm{IMPSS}}$ and the optimal IMPSS preconditioner ${\mathcal{P}}_{\mathrm{OIMPSS}}$ be defined as in equations (7) and (10), respectively. Then there are $n+i$ $(0\le i\le m)$ linearly independent eigenvectors of the preconditioned matrix ${\mathcal{P}}_{\mathrm{IMPSS}}^{-1}\mathcal{A}$ :

• n eigenvectors

$$\left[\begin{array}{c}{u}_p^1\hfill \\ \multicolumn{1}{c}{0}\end{array}\right]$$

( $p=1,2,\dots ,n$ ) that correspond to the eigenvalue $1$ , where $u_p^1$ $(p=1,2,\dots ,n)$ are arbitrary linearly independent vectors;

• i ( $0\le i\le m$ ) eigenvectors

$$\left[\begin{array}{c}{u}_p^2\hfill \\ \multicolumn{1}{c}{v_p^2}\end{array}\right]$$

( $0\le p\le i$ ) that correspond to nonunit eigenvalues, where $v_p^2(\ne 0)$ satisfies $(C+B{A}^{-1}{B}^{*})v_p^2=\lambda \hat{S}{v}_p^2$ and $u_p^2=[(\lambda /\alpha (1-\lambda ))I-A^{-1}]B^{*}{v}_p^2$ .

In addition, the preconditioned matrix ${\mathcal{P}}_{\mathrm{OIMPSS}}^{-1}\mathcal{A}$ has only n linearly independent eigenvectors

$$\left[\begin{array}{c}{u}_p^1\hfill \\ \multicolumn{1}{c}{0}\end{array}\right]$$

( $p=1,2,\dots ,n$ ), where $u_p^1$ ( $p=1,2,\dots ,n$ ) are arbitrary linearly independent vectors.

Proof

We first study the eigenvector distribution of the IMPSS preconditioned matrix ${\mathcal{P}}_{\mathrm{IMPSS}}^{-1}\mathcal{A}$ . Let $\lambda$ be an eigenvalue of ${\mathcal{P}}_{\mathrm{IMPSS}}^{-1}\mathcal{A}$ and

$$\left[\begin{array}{c}u\hfill \\ \multicolumn{1}{c}{v}\end{array}\right]$$

be the corresponding eigenvector. Then from equation (8), we have the following eigenvalue problem

$$\begin{array}{l}\left[\begin{array}{cc}I& A^{-1}{B}^{*}-\left(A^{-1}+\frac{1}{\alpha }I\right)B^{*}{\widehat{S}}^{-1}(C+B{A}^{-1}{B}^{*})\\ 0& {\widehat{S}}^{-1}(C+B{A}^{-1}{B}^{*})\end{array}\right]\\ \left[\begin{array}{c}u\\ v\end{array}\right]=\lambda \left[\begin{array}{c}u\\ v\end{array}\right]\end{array}$$

(11)

Expanding equation (11) out, we have

$$(1-\lambda )u+(1-\lambda )A^{-1}{B}^{*}v-\frac{\lambda }{\alpha }{B}^{*}v=0$$

(12)

and

$$(C+B{A}^{-1}{B}^{*})v=\lambda \hat{S}v$$

(13)

From the eigenvalue distribution observed in Theorem 1, we now consider two cases. If $\lambda =1$ , then equation (12) becomes $B^{*}v=0$ , which results in $v=0$ since $B^{*}$ is of full column rank. Let $u_p^1\in {\mathbb{C}}^n$ ( $p=1,\dots ,n$ ) be arbitrary linearly independent vectors. Then there are n eigenvectors of the form

$$\left[\begin{array}{c}{u}_p^1\hfill \\ \multicolumn{1}{c}{0}\end{array}\right]$$

that correspond to the eigenvalue $1$ .

If $\lambda \ne 1$ , then we have from equation (12)

$$u=[\frac{\lambda }{\alpha (1-\lambda )}I-A^{-1}]B^{*}v$$

(14)

and the vector v satisfies equation (13). Although equation (13) is trivially satisfied by $v=0$ , this cannot happen. Otherwise, from equation (14), we have $u=0$ , a contradiction. Therefore, if there exists any $v\ne 0$ which satisfies equation (13), then there are i ( $0\le i\le m$ ) eigenvectors

$$\left[\begin{array}{c}{u}_p^2\hfill \\ \multicolumn{1}{c}{v_p^2}\end{array}\right]$$

( $0\le p\le i$ ) that correspond to nonunit eigenvalues, where $v_p^2(\ne 0)$ satisfies $(C+B{A}^{-1}{B}^{*})v_p^2=\lambda \hat{S}{v}_p^2$ and $u_p^2=[(\lambda /\alpha (1-\lambda ))I-A^{-1}]B^{*}{v}_p^2$ .

Now, we consider the eigenvector distribution of the preconditioned matrix ${\mathcal{P}}_{\mathrm{OIMPSS}}^{-1}\mathcal{A}$ . From equation (9), we have

$$\left[\begin{array}{cc}I\hfill & -\frac{1}{\alpha }{B}^{*}\hfill \\ \multicolumn{1}{c}{0}& I\hfill \end{array}\right]\left[\begin{array}{c}u\hfill \\ \multicolumn{1}{c}{v}\end{array}\right]=\left[\begin{array}{c}u\hfill \\ \multicolumn{1}{c}{v}\end{array}\right]$$

which is equivalent to the following two equalities

$$u-\frac{1}{\alpha }{B}^{*}v=u\mathrm{and}v=v$$

The second equality is satisfied naturally. From the first equality, we know that $v=0$ since the matrix $B^{*}$ has full column rank. Therefore, the OIMPSS preconditioned matrix ${\mathcal{P}}_{\mathrm{OIMPSS}}^{-1}\mathcal{A}$ has only n linearly independent eigenvectors

$$\left[\begin{array}{c}{u}_p^1\hfill \\ \multicolumn{1}{c}{0}\end{array}\right]$$

( $p=1,2,\dots ,n$ ), where $u_p^1\in {\mathbb{C}}^n$ ( $p=1,2,\dots ,n$ ) are arbitrary linearly independent vectors.

The linear independence of these eigenvectors can be derived by the method studied in Cao et al. ¹⁰ [Theorem 3.2] with only technical modifications. Hence, we omit the rest of the proof here.

Based on the special structures of the preconditioned matrices ${\mathcal{P}}_{\mathrm{IMPSS}}^{-1}\mathcal{A}$ and ${\mathcal{P}}_{\mathrm{OIMPSS}}^{-1}\mathcal{A}$ , the properties with respect to the minimal polynomials of these preconditioned matrices are given in Theorem 3.

Theorem 3

Let $A\in {\mathbb{C}}^{n\times n}$ be a non-Hermitian positive-definite matrix, $B\in {\mathbb{R}}^{m\times n}$ have full row rank, $C\in {\mathbb{R}}^{m\times m}$ be a Hermitian positive-semidefinite matrix. Let $\hat{S}\in {\mathbb{R}}^{m\times m}$ be a nonsingular matrix and $\alpha$ be a positive constant. Let the IMPSS preconditioner ${\mathcal{P}}_{\mathrm{IMPSS}}$ and the optimal IMPSS preconditioner ${\mathcal{P}}_{\mathrm{OIMPSS}}$ be defined as in equations (7) and (10), respectively. Then, both the degree of the minimal polynomial of the preconditioned matrix ${\mathcal{P}}_{\mathrm{IMPSS}}^{-1}\mathcal{A}$ and the dimension of the Krylov subspace $\mathcal{K}({\mathcal{P}}_{\mathrm{IMPSS}}^{-1}\mathcal{A},b)$ are at most $m+1$ . In addition, the minimal polynomial of the preconditioned matrix ${\mathcal{P}}_{\mathrm{OIMPSS}}^{-1}\mathcal{A}$ is

$$({\mathcal{P}}_{\mathrm{OIMPSS}}^{-1}\mathcal{A}-I{)}^2$$

(15)

Proof

Let ${\mathrm{\Delta }}_1=A^{-1}{B}^{\ast }-(A^{-1}+(1/\alpha )I)B^{\ast }{\widehat{S}}^{-1}(C+B{A}^{-1}{B}^{\ast })$ and ${\mathrm{\Delta }}_2={\hat{S}}^{-1}(C+B{A}^{-1}{B}^{*})$ . Then the preconditioned matrix ${\mathcal{P}}_{\mathrm{IMPSS}}^{-1}\mathcal{A}$ can be rewritten as the block upper triangular matrix

$${\mathcal{P}}_{\mathrm{IMPSS}}^{-1}\mathcal{A}=\left[\begin{array}{cc}I\hfill & {\mathrm{\Delta }}_1\hfill \\ \multicolumn{1}{c}{0}& {\mathrm{\Delta }}_2\hfill \end{array}\right]$$

Let ${\lambda }_i$ ( $i=1,\dots ,m$ ) be the eigenvalues of the $m\times m$ matrix ${\mathrm{\Delta }}_2$ . Then ${\lambda }_i$ are also the eigenvalues of the preconditioned matrix ${\mathcal{P}}_{\mathrm{IMPSS}}^{-1}\mathcal{A}$ . According to the eigenvalue distribution studied in Theorem 1, the characteristic polynomial of the preconditioned matrix ${\mathcal{P}}_{\mathrm{IMPSS}}^{-1}\mathcal{A}$ can be defined as

$${\mathrm{\Phi }}_{{\mathcal{P}}_{IMPSS}^{-1}\mathcal{A}}(\lambda )=det({\mathcal{P}}_{IMPSS}^{-1}\mathcal{A}-\lambda I)=(-{1)}^{m+n}{(\lambda -1)}^n{\displaystyle \prod _{i=1}^m}(\lambda -{\lambda }_i)$$

where $\det (·)$ denotes the determinant of the corresponding matrix. Consider the following polynomial of degree $m+1$

$$\mathrm{\Psi }(\lambda )=(\lambda -1)\underset{i=1}{\overset{m}{\Pi }}(\lambda -{\lambda }_i)$$

Then we have

$$\begin{array}{l}\mathrm{\Psi }({\mathcal{P}}_{\mathcal{I}\mathcal{M}\mathcal{P}\mathcal{S}\mathcal{S}}^{-\infty }\mathcal{A})=({\mathcal{P}}_{IMPSS}^{-1}\mathcal{A}-I){\displaystyle \prod _{i=1}^m({\mathcal{P}}_{IMPSS}^{-1}\mathcal{A}-{\lambda }_{i}I)}\\ =\left[\begin{array}{cc}0& {\mathrm{\Delta }}_1{\displaystyle \prod _{i=1}^m({\mathrm{\Delta }}_2-{\lambda }_{i}I)}\\ 0& ({\mathrm{\Delta }}_2-I){\displaystyle \prod _{i=1}^m({\mathrm{\Delta }}_2-{\lambda }_{i}I)}\end{array}\right]\end{array}$$

By the Hamilton–Cayley theorem, we have $\underset{i=1}{\overset{m}{\Pi }}({\mathrm{\Delta }}_2-{\lambda }_{i}I)=0$ . Thus, the degree of the minimal polynomial of the preconditioned matrix ${\mathcal{P}}_{\mathrm{IMPSS}}^{-1}\mathcal{A}$ is at most $m+1$ . From Saad ² [Proposition 6.1], we know that the dimension of the Krylov subspace $\mathcal{K}({\mathcal{P}}_{\mathrm{IMPSS}}^{-1}\mathcal{A},b)$ is at most $m+1$ .

From equation (9), we can see that the preconditioned matrix ${\mathcal{P}}_{\mathrm{OIMPSS}}^{-1}\mathcal{A}$ is a block upper triangular matrix with block diagonal matrices being the identity matrices. Therefore, the corresponding minimal polynomial (15) of the preconditioned matrix ${\mathcal{P}}_{\mathrm{OIMPSS}}^{-1}\mathcal{A}$ can be obtained.

Numerical experiments

In this section, we present some numerical examples to illustrate the effectiveness of the new proposed IMPSS preconditioners. To this end, we will use the IMPSS preconditioned GMRES iteration method to solve the generalized saddle point problems (1) arising from the $Q1-P0$ mixed finite element discretization of the following two-dimensional steady incompressible linearized Navier–Stokes equation ²⁷

$$\{\begin{array}{c}-\nu \mathrm{\Delta }u+\omega ·\nabla u+\nabla p=f\hfill \\ \multicolumn{1}{c}{\nabla ·u=0}\end{array}\mathrm{in}\mathrm{\Omega }$$

(16)

where $\mathrm{\Omega }=[0,1]\times [0,1]$ is the unit bounded square domain, $\nu >0$ is the viscosity, $\omega$ is the velocity field obtained from the previous Picard step, u and p stand for the velocity and the pressure, respectively. The test model is the lid-driven cavity problem. The IFISS software package written by Elman et al. ²⁷ is used to generate the test generalized saddle point problems. Here, we choose uniform grids of increasing sizes ( $8\times 8$ , $16\times 16$ , $32\times 32$ , and $64\times 64$ ) to discretize the problem domain and take the local stabilization method for the $Q1-P0$ element (the stabilization parameter is taken as 0.25).

To emphasize the distinct advantages of these new IMPSS preconditioners, the GMRES iteration methods with the original DPSS preconditioner (3) and its improved variant MPSS preconditioner (4) are used to solve test problems. In particular, the left preconditioning technique is used. We compare these preconditioners by the number of iteration steps (denoted by “IT”) and the elapsed CPU times in seconds (denoted by “CPU”). In actual computations, the initial guess is the zero vector and the stopping criterion is a reduction of at least six orders of magnitude of the initial residual norm (denoted by “RES”). At each iteration of the IMPSS, the DPSS, the MPSS preconditioned GMRES method, some sub-systems of linear equations need to be solved. If the coefficient matrix is non-Hermitian positive-definite, then we use the sparse LU factorization method. If the coefficient matrix is Hermitian positive-definite, then we choose the sparse Cholesky factorization method. All trial runs are performed in MATLAB2010 on an Intel Core (4G RAM) Windows 7 system.

As discussed in section “Spectral analysis of the IMPSS preconditioned matrix,” the preconditioning matrix $\hat{S}$ plays an important role in the IMPSS preconditioners. We need to choose the preconditioning matrix $\hat{S}$ with cheapest computational cost and best convergence performance. It is shown in Wathen and Silvester ²⁸ and Wathen and Silvester ²⁹ that for the Stokes equations, there exists a pressure mass matrix $M_p$ (i.e. the Grammian matrix of basis functions for discrete pressure subspace) such that

$${\gamma }^2\le \frac{v^{*}(C+B{A}^{-1}{B}^{*})v}{v^{*}{M}_{p}v}\le {\mathrm{\Gamma }}^2,\forall v\in {\mathbb{C}}^m\backslash \{0\}$$

where constants $\gamma$ and $\mathrm{\Gamma }$ are independent of the meshsize h. For the $Q1-P0$ stabilized finite element method, the diagonal matrix $D_c=h^{d}I$ with d be the spatial dimension is spectral equivalent to the pressure mass matrix, that is, there exist constants $\theta$ , $\mathrm{\Theta }$ independent of the meshsize h such that

$${\theta }^2\le \frac{v^{*}{M}_{p}v}{v^{*}{D}_{c}v}\le {\mathrm{\Theta }}^2,\forall v\in {\mathbb{C}}^m\backslash \{0\}$$

Hence, we can choose $\hat{S}=D_c$ for the IMPSS preconditioner such that ${\hat{S}}^{-1}$ is very easy to obtain. The theoretical optimal OIMPSS preconditioner ( $\hat{S}=C+B{A}^{-1}{B}^{*}$ ) is also studied. In addition, the parameters used in the discussed preconditioners should be chosen appropriately. In general, finding these optimal parameters is a very difficult task. Here, we use some algebraic estimation methods studied in Cao et al. ¹⁴ and Ren and Cao, ³⁰ which can also lead to good numerical results. Specifically, the parameter estimate method proposed in Ren and Cao ³⁰ is used to derive the following estimation for the parameter of the DPSS preconditioner

$${\alpha }_{\mathrm{DPSS}}=\frac{\parallel \mathcal{P}{\parallel }_F+\parallel \mathcal{S}{\parallel }_F}{2(n+m)}$$

(17)

which is the minimum value of the following function with respect to the parameter $\alpha$

$$\mathrm{\Phi }(\alpha )=(n+m){\alpha }^2-\alpha (\parallel \mathcal{P}{\parallel }_F+\parallel \mathcal{S}{\parallel }_F)+\parallel \mathcal{P}{\parallel }_F\parallel \mathcal{S}{\parallel }_F$$

The basic idea of constructing the above quadratic equation is to obtain a practical parameter ${\alpha }_{\mathrm{DPSS}}$ such that the DPSS preconditioner approximates the generalized saddle point matrix as mush as possible in certain norm. For details of the algebraic estimate method, see Ren and Cao ³⁰ [Section 4]. Here, the matrices $\mathcal{P}$ and $\mathcal{S}$ are defined as in equation (2). We choose ${\alpha }_{\mathrm{MPSS}}=\parallel B{B}^{*}{\parallel }_2/\parallel C{\parallel }_2$ (which balances the matrices B and C) ¹⁴ and ${\beta }_{\mathrm{MPSS}}=0.25$ (which is the same as the local stabilization parameter) for the MPSS preconditioner. Here, $\parallel ·{\parallel }_F$ and $\parallel ·{\parallel }_2$ denote the Frobenius norm and the Euclidean norm of the corresponding matrix, respectively. According to the eigenvalue distribution of the IMPSS preconditioned matrix (see Theorem 1) and the minimal polynomial of the theoretical optimal IMPSS preconditioned matrix (see Theorem 3), we can always take $\alpha =1.0$ for the IMPSS preconditioner and the OIMPSS preconditioner since the parameter $\alpha$ may not greatly impact the preconditioning effects. ¹⁷

In Table 1, we list the estimated parameters as well as the numerical results of the preconditioned GMRES iteration methods for $\nu =1$ . In Table 2, we list the same items for $\nu =0.1$ . In these tables, numerical results of the standard GMRES method (denoted by I) are listed to show the advantages of the preconditioning effects. In addition, to show the clustering property of the eigenvalues of the preconditioned matrices, using the respective estimated parameters listed in Tables 1 and 2, we plot the eigenvalue distribution of the generalized saddle point matrices, the DPSS preconditioned matrices, the MPSS preconditioned matrices, and the IMPSS preconditioned matrices for the $16\times 16$ mesh grids with $\nu =1$ and $\nu =0.1$ in Figures 1 and 2, respectively. In these figures, the horizontal ordinate and the longitudinal coordinate indicate the real part and the imaginary part of the eigenvalues, respectively.

Table 1.

Numerical results of the preconditioned GMRES method for Oseen problem with $\nu =1$ .

Prec.		$8\times 8$	$16\times 16$	$32\times 32$	$64\times 64$
I	IT	77	171	368	404
	CPU	0.0635	0.9208	8.3449	32.5204
	RES	7.9934e–7	9.8576e–7	9.4717e–7	9.9803e–7
DPSS	$\alpha$	0.0721	0.0388	0.0201	0.0102
	IT	55	70	86	115
	CPU	0.0242	0.0679	0.3196	3.1089
	RES	7.8890e–7	8.7321e–7	8.6436e–7	9.5772e–7
MPSS	$(\alpha ,\beta )$	(3.7038, 0.25)	(3.9237, 0.25)	(3.9685, 0.25)	(3.9670, 0.25)
	IT	17	22	25	29
	CPU	0.0073	0.2040	3.5245	183.6651
	RES	2.7112e–7	7.0980e–7	4.4526e–7	9.3085e–7
OIMPSS	IT	3	3	3	3
	CPU	0.0034	0.0603	2.3664	172.1122
	RES	1.1968e–15	9.1926e–15	6.5605e–14	2.5373e–13
IMPSS	IT	24	28	27	22
	CPU	0.0081	0.0329	0.2325	0.6001
	RES	7.6281e–7	4.5571e–7	4.5911e–7	7.0760e–7

DPSS: deteriorated positive-definite and skew-Hermitian splitting; GMRES: generalized minimal residual; MPSS: modified positive-definite and skew-Hermitian splitting; OIMPSS: optimal modified positive-definite and skew-Hermitian splitting; IMPSS: inexact modified positive-definite and skew-Hermitian splitting.

Table 2.

Numerical results of the preconditioned GMRES method for Oseen problem with $\nu =0.1$ .

Prec.		$8\times 8$	$16\times 16$	$32\times 32$	$64\times 64$
I	IT	84	157	275	539
	CPU	0.0794	0.3316	3.0207	28.9593
	RES	8.9902e–7	9.0207e–7	9.5096e–7	9.7511e–7
DPSS	$\alpha$	0.0268	0.0099	0.0037	0.0015
	IT	36	56	81	114
	CPU	0.0155	0.1449	0.4603	3.0747
	RES	9.1167e–7	6.9278e–7	8.2637e–7	9.4520e–7
MPSS	$(\alpha ,\beta )$	(3.7038, 0.25)	(3.9237, 0.25)	(3.9685, 0.25)	(3.9670, 0.25)
	IT	10	15	18	20
	CPU	0.0054	0.0732	3.4272	182.4161
	RES	3.4673e–7	2.6273e–7	5.4827e–7	3.6449e–7
OIMPSS	IT	3	3	3	3
	CPU	0.0045	0.0634	2.2084	172.3443
	RES	2.5557e–15	1.3322e–14	6.4613e–14	2.1516e–13
IMPSS	IT	18	19	19	19
	CPU	0.0064	0.0198	0.1841	0.6887
	RES	4.9360e–7	4.2812e–7	8.9734e–7	9.3484e–7

DPSS: deteriorated positive-definite and skew-Hermitian splitting; MPSS: modified positive-definite and skew-Hermitian splitting; OIMPSS: optimal modified positive-definite and skew-Hermitian splitting; IMPSS: inexact modified positive-definite and skew-Hermitian splitting.

Figure 1.

Eigenvalue distribution of the preconditioned matrices ( $16\times 16$ grids, $\nu =1$ ).

Figure 2.

Eigenvalue distribution of the preconditioned matrices ( $16\times 16$ grids, $\nu =0.1$ ).

From the numerical results in Tables 1 and 2, we can see that without preconditioning, the GMRES method converges very slow. Just considering iteration counts, the discussed preconditioners can greatly accelerate the convergence rate of the GMRES method. However, the elapsed CPU times show that not all discussed preconditioners are very efficient, especially for large problems. Although the iteration counts of the MPSS preconditioned methods are less than those of the DPSS preconditioned methods, the elapsed CPU times show that the MPSS preconditioner requires a much higher cost, in terms of iteration counts, than the original DPSS preconditioner. So, the MPSS preconditioner does not improve the DPSS preconditioner. Using inexact preconditioning techniques, our new IMPSS preconditioner not only improves the DPSS preconditioner but also dramatically improves the computing efficiency of the MPSS preconditioner. In addition, from these numerical results, we can see that for the test problems given in equation (16), the iteration counts of the IMPSS preconditioned method do not change dramatically for different meshsizes when the viscosity parameter $\nu =1$ , the IMPSS preconditioner even shows an h-independent convergence property when the viscosity parameter $\nu =0.1$ . The OIMPSS preconditioned GMRES iteration methods converge to the exact solution within three iteration steps for all test problems. Therefore, from the iteration counts of these preconditioned GMRES methods, the OIMPSS preconditioner is superior and shows an h-independent convergence property, which confirms our theoretical results in Theorems 2 and 3. However, the OIMPSS preconditioner is also very computationally costly. Figures 1 and 2 show that all preconditioned matrices have clustered eigenvalue distribution, which is in agreement with the theoretical results in Theorem 1. Therefore, we can conclude that the new proposed IMPSS preconditioner is very efficient and greatly improves the DPSS and the MPSS preconditioners.