A block-diagonally preconditioned Uzawa splitting iteration method for solving a class of saddle-point problems

Abstract

This paper develops a block diagonal preconditioned Uzawa splitting (BDP-US) method for solving saddle point problems by generalizing the Uzawa splitting iteration method proposed by Li and Ma (Numer Math Theory Methods Appl 2018; 11: 235–246). A sufficient condition is then provided to ensure the convergence of the BDP-US method. Meanwhile, a preconditioner on the basis of the BDP-US method is proposed, the spectral properties of the preconditioned matrix is analyzed, and the choice of the parameters for this matrix splitting iteration method is discussed. Numerical results are provided to support the obtained results, and demonstrate the effectiveness of BDP-US method as well as the corresponding preconditioner.

Mathematics Subject Classification: 65F10 65F08

Introduction

Consider the linear system of the $2\times 2$ block structure:

$$\mathcal{A}z\equiv \left(\begin{array}{cc}A\hfill & B\hfill \\ \multicolumn{1}{c}{-B^T}& O\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{R}}^{m\times m}$ is symmetric and positive definite, $B\in {\mathbb{R}}^{m\times n}$ with $n<<m$ , $x,f\in {\mathbb{R}}^m$ , $y,g\in {\mathbb{R}}^n$ , $m=2p^2,n=p^2$ , and $p$ denotes the mesh size. $B^T$ is the transpose of $B$ , and $O$ is the zero matrix of proper dimension. It is known that the linear system has a unique solution when $B$ is of full column rank (see Benzi et al., ¹ Bai and Bai, ² Bai and Pan ³ ).

Linear system of form (1) is called a saddle-point problem, and arises from many important applications in various fields of science and engineering, including computational fluid dynamics, optimization, optimal control, constrained and weighted least-squares estimation as well as many others.^4–8 Then, designing effective algorithms to solve it has attracted researchers’ attention. In general, for a large scale problem (1), only iterative methods are computationally feasible due to excessive storage and/or computational cost. By using some algebraic properties and the structure of the coefficient matrix, many iterative methods have been developed to solve (1), such as Uzawa type methods,^9–12 SOR (successive overrelaxation) type methods,^13–16 HSS (Hermitian and skew-Hermitian splitting) type methods and its accelerated variants,^17–20,28 and Krylov subspace methods with suitable preconditioners.^21–24 In particular, Bai ²⁸ discussed the optimal parameters in the HSS-like methods for saddle-point problems, and provided the quasi-optimal parameters of the HSS iteration method.

Recently, based on the shift-splitting iteration method ²¹ and the classical Uzawa method, ⁹ Li and Ma ²⁵ proposed the following Uzawa splitting (US) iteration method

$$\{\begin{array}{c}{x}^{k+1}=x^k+2(\alpha I_m+A{)}^{-1}(f-A{x}^k-B{y}^k)\hfill \\ \multicolumn{1}{c}{y^{k+1}=y^k+\tau Q^{-1}(B^T{x}^{k+1}-g)}\end{array}$$

to solve singular saddle-point problem (1) with $\mathrm{r}\mathrm{ank}(B)\le n$ and $\mathrm{r}\mathrm{ank}(·)$ being the rank of the matrix, where $Q\in {\mathbb{R}}^{n\times n}$ is a symmetric positive definite matrix, $\alpha$ and $\tau$ are two positive constants. Even though this method is a special case of the parameterized inexact Uzawa (PIU) method, ¹⁵ it has better performance than the Uzawa method and Uzawa-SOR methods. ²⁶ Moreover, this method is convergent when

$$\frac{2\alpha }{\tau }>{\lambda }_{max},$$

where ${\lambda }_{max}$ is the largest eigenvalue of the matrix $B{Q}^{-1}{B}^T$ . But, the US method might converge very slowly as the scale of the problem increases. In fact, at this case, the value of ${\lambda }_{max}$ may become larger, which makes $\alpha$ too large, or $\tau$ too small. For example, when $p=8,16,24$ , the given values of $(\alpha ,\tau )$ in Li and Ma ²⁵ are $(221,0.57)$ , $(921,0.47)$ , and $(1237,0.46)$ , respectively. It should be mentioned that no strategy for selecting iteration parameters can be found in Li and Ma, ²⁵ and it requires to accurately solve a linear system with the coefficient matrix $\alpha I_m+A$ in per iteration of the US method. However, it is usually impractical to directly solve (1), and may be expensive to solve it by the conjugate gradient method without preconditioning when $m$ is large.

To overcome the above shortfalls, in this paper, we (a) develop a BDP-US method for solving the saddle-point problem (1) by generalizing the US method; (b) provide a sufficient condition to ensure the convergence of the BDP-US method; (c) present a splitting preconditioner and analyze the spectral properties of the preconditioned matrix; and (d) discuss the choice of the parameters for this BDP-US method. The proposed method includes the US method in Li and Ma ²⁵ as a special case. Numerical results demonstrate that the BDP-US method is effective and feasible for solving the singular saddle point problems (1). Moreover, the splitting preconditioner can improve the convergence rate of Generalized Minimal Residual (GMRES) method.

The remainder of this paper is organized as follows. In Section 2, the BDP-US iteration method is presented. The convergence properties of the BDP-US method and the spectral properties of the preconditioned matrix are analyzed in Section 3 and 4, respectively. The choice of the parameters for the BDP-US method is discussed in Section 5, and numerical experiments are provided in Section 6. Finally, we draw some conclusions in Section 7.

In the following, we denote by $I_n$ the identity matrix of order $n$ , by $\mathrm{Re}(·)$ the real part of a complex number, by the superscripts $T$ the transpose of a matrix or vector, by $\rho (·)$ the spectral radius of a matrix, by $(u;v)$ the vector $(u^{*},v^{*}{)}^{*}$ and by ${\lambda }_{max}(·)$ the largest eigenvalues of a symmetric and positive definite matrix.

The BDP-US iteration method

In this section, we shall develop the BDP-US iteration method for solving the linear system (1).

Let

$$\begin{array}{c}{P}^{-1}=\left(\begin{array}{cc}{P}_1\hfill & O\hfill \\ \multicolumn{1}{c}{O}& P_2\hfill \end{array}\right),\hfill \end{array}$$

where $P_1\in {\mathbb{R}}^{m\times m}$ and $P_2\in {\mathbb{R}}^{n\times n}$ are positive definite matrices. Then by left and right multiplying the linear system (1) with $P^{-1}$ and $P^{-T}$ respectively, we obtain

$$\begin{array}{c}\left(\begin{array}{cc}{P}_1\hfill & O\hfill \\ \multicolumn{1}{c}{O}& P_2\hfill \end{array}\right)\left(\begin{array}{cc}A\hfill & B\hfill \\ \multicolumn{1}{c}{-B^T}& O\hfill \end{array}\right)\left(\begin{array}{cc}{P}_1^T\hfill & O\hfill \\ \multicolumn{1}{c}{O}& P_2^T\hfill \end{array}\right)\left(\begin{array}{cc}{P}_1^{-T}\hfill & O\hfill \\ \multicolumn{1}{c}{O}& P_2^{-T}\hfill \end{array}\right)\left(\begin{array}{c}x\hfill \\ \multicolumn{1}{c}{y}\end{array}\right)\hfill \\ =\left(\begin{array}{cc}{P}_1\hfill & O\hfill \\ \multicolumn{1}{c}{O}& P_2\hfill \end{array}\right)\left(\begin{array}{c}f\hfill \\ \multicolumn{1}{c}{-g}\end{array}\right),\hfill \end{array}$$

which can be written as

$$\hat{\mathcal{A}}\hat{z}\equiv \left(\begin{array}{cc}{P}_1{\mathrm{AP}}_1^T\hfill & P_1{\mathrm{BP}}_2^T\hfill \\ \multicolumn{1}{c}{-P_2{B}^T{P}_1^T}& O\hfill \end{array}\right)\left(\begin{array}{c}{P}_1^{-T}x\hfill \\ \multicolumn{1}{c}{P_2^{-T}y}\end{array}\right)=\left(\begin{array}{c}{P}_{1}f\hfill \\ \multicolumn{1}{c}{-P_{2}g}\end{array}\right)\equiv \hat{b}.$$

(2)

Similar to Bai et al., ²¹ we can define the corresponding shift-splitting of $P_1{\mathrm{AP}}_1^T$ as

$$\begin{array}{c}{P}_1{\mathrm{AP}}_1^T=\frac{1}{2}(\alpha I_m+P_1{\mathrm{AP}}_1^T)-\frac{1}{2}(\alpha I_m-P_1{\mathrm{AP}}_1^T),\hfill \end{array}$$

(3)

where $\alpha$ is a positive constant. For the available approximation ${\hat{z}}^k=((P_1^{-T}{x}^k{)}^T,(P_2^{-T}{y}^k{)}^T{)}^T$ , $k=0,1,\dots$ , the following iteration scheme is first applied to the first $m$ equations in (2) to compute the update $x^{k+1}$ :

$$\begin{array}{c}\multicolumn{1}{c}{\frac{1}{2}(\alpha I_m+P_1{\mathrm{AP}}_1^T)(P_1^{-T}{x}^{k+1})=\frac{1}{2}(\alpha I_m-P_1{\mathrm{AP}}_1^T)(P_1^{-T}{x}^k)}\\ +(P_{1}f-P_1{\mathrm{BP}}_2^T(P_2^{-T}{y}^k)).\hfill \end{array}$$

(4)

Then the Uzawa iteration method is used to the last $n$ equations in (2) to compute the update $y^{k+1}$ :

$$\frac{1}{\tau }Q(P_2^{-T}{y}^{k+1}-P_2^{-T}{y}^k)=P_2{B}^T{P}_1^T(P_1^{-T}{x}^{k+1})-P_{2}g,$$

(5)

where $Q\in {\mathbb{R}}^{n\times n}$ is symmetric and positive definite, and $\tau$ is a positive constant.

For simplicity, we denote $\hat{x}=P_1^{-T}x$ , $\hat{y}=P_2^{-T}y$ , ${\hat{x}}^k=P_1^{-T}{x}^k$ and ${\hat{y}}^k=P_2^{-T}{y}^k$ . Then the results in (4) and (5) suggest the following BDP-US method for the linear system (2):

$$\begin{array}{c}\multicolumn{1}{c}{\{\begin{array}{c}{\hat{x}}^{k+1}={\hat{x}}^k+2(\alpha I_m+P_1{\mathrm{AP}}_1^T{)}^{-1}(P_{1}f-P_1{\mathrm{AP}}_1^T{\hat{x}}^k-P_1{\mathrm{BP}}_2^T{\hat{y}}^k),\hfill \\ \multicolumn{1}{c}{{\hat{y}}^{k+1}={\hat{y}}^k+\tau Q^{-1}(P_2{B}^T{P}_1^T{\hat{x}}^{k+1}-P_{2}g),}\end{array}}\end{array}$$

(6)

or

$$\begin{array}{c}{\hat{z}}^{k+1}=T(\alpha ,\tau ){\hat{z}}^k+H(\alpha ,\tau )\hat{b},\hfill \end{array}$$

(7)

where the iteration matrix

$$T(\alpha ,\tau )=\left(\begin{array}{cc}{T}_{11}(\alpha ,\tau )\hfill & T_{12}(\alpha ,\tau )\hfill \\ \multicolumn{1}{c}{T_{21}(\alpha ,\tau )}& T_{22}(\alpha ,\tau )\hfill \end{array}\right)$$

with

$$\begin{array}{c}{T}_{11}(\alpha ,\tau )=(\alpha I_m+P_1{\mathrm{AP}}_1^T{)}^{-1}(\alpha I_m-P_1{\mathrm{AP}}_1^T),\hfill \\ T_{12}(\alpha ,\tau )=-2(\alpha I_m+P_1{\mathrm{AP}}_1^T{)}^{-1}(P_1{\mathrm{BP}}_2^T),\hfill \\ T_{21}(\alpha ,\tau )=\tau Q^{-1}(P_2{B}^T{P}_1^T)(\alpha I_m+P_1{\mathrm{AP}}_1^T{)}^{-1}(\alpha I_m-P_1{\mathrm{AP}}_1^T),\hfill \\ T_{22}(\alpha ,\tau )=I_n-2\tau Q^{-1}(P_2{B}^T{P}_1^T)(\alpha I_m+P_1{\mathrm{AP}}_1^T{)}^{-1}(P_1{\mathrm{BP}}_2^T),\hfill \end{array}$$

and

$$H(\alpha ,\tau )=\left(\begin{array}{cc}2(\alpha I_m+P_1{\mathrm{AP}}_1^T{)}^{-1}\hfill & O\hfill \\ \multicolumn{1}{c}{2\tau Q^{-1}(P_2{B}^T{P}_1^T)(\alpha I_m+P_1{\mathrm{AP}}_1^T{)}^{-1}}& \tau Q^{-1}\hfill \end{array}\right).$$

On the other hand, let

$$M(\alpha ,\tau )=\left(\begin{array}{cc}\frac{1}{2}(\alpha I_m+P_1{\mathrm{AP}}_1^T)\hfill & O\hfill \\ \multicolumn{1}{c}{-P_2{B}^T{P}_1^T}& \frac{1}{\tau }Q\hfill \end{array}\right)$$

(8)

and

$$N(\alpha ,\tau )=\left(\begin{array}{cc}\frac{1}{2}(\alpha I_m-P_1{\mathrm{AP}}_1^T)\hfill & -P_1{\mathrm{BP}}_2^T\hfill \\ \multicolumn{1}{c}{O}& \frac{1}{\tau }Q\hfill \end{array}\right),$$

(9)

then we obtain a splitting $\mathcal{A}=M(\alpha ,\tau )-N(\alpha ,\tau )$ , and

$$T(\alpha ,\tau )=M(\alpha ,\tau {)}^{-1}N(\alpha ,\tau ).$$

(10)

Overall, we can perform the BDP-US method as follows: Given an initial guess ${\hat{z}}^0=(({\hat{x}}^0{)}^T,({\hat{y}}^0{)}^T{)}^T$ with ${\hat{x}}^0\in {\mathbb{R}}^m$ and ${\hat{y}}^0\in {\mathbb{R}}^n$ , two positive parameter $\alpha$ and $\tau$ ; for $k=0,1,2,\dots$ , compute ${\hat{z}}^{k+1}$ by the iteration scheme (6) until the iteration sequence $\{{\hat{z}}^k{\}}_{k=0}^{\infty }$ converges. Thus we obtain the converged sequence $\{(x^k,y^k)\}$ with $x^k=P_1^T{\hat{x}}^k$ and $y^k=P_1^T{\hat{y}}^k$ .

Clearly, this BDP-US method reduces to the US method in Li and Ma ²⁵ by setting $P_1=I_m$ and $P_2=I_n$ , and different $Q$ , $\alpha$ and $\tau$ can be selected in actual computation.

Convergence analysis of the BDP-US method

In this section, we shall analyze the convergence of the BDP-US iteration method, and present a sufficient condition for its convergence. To this end, we need the following result.

• Lemma 3.1 ²⁷ Both roots of the real quadratic equation $x^2-\mathrm{bx}+c=0$ are less than one in modulus if and only if $|c|<1$ and $|b|<1+c$ .

It is well known that the BDP-US method (7) is convergent if and only if $\rho (T(\alpha ,\tau ))<1$ . Let $\lambda$ be an eigenvalue of $T(\alpha ,\tau )$ , and $(u;v)$ with $u\in {\mathbb{C}}^m$ and $v\in {\mathbb{C}}^n$ be the corresponding eigenvector, then

$$T(\alpha ,\tau )\left(\begin{array}{c}u\hfill \\ \multicolumn{1}{c}{v}\end{array}\right)=\lambda \left(\begin{array}{c}u\hfill \\ \multicolumn{1}{c}{v}\end{array}\right).$$

Namely,

$$N(\alpha ,\tau )\left(\begin{array}{c}u\hfill \\ \multicolumn{1}{c}{v}\end{array}\right)=\lambda M(\alpha ,\tau )\left(\begin{array}{c}u\hfill \\ \multicolumn{1}{c}{v}\end{array}\right).$$

This, (8) and (9) imply that

$$\begin{array}{c}\left(\begin{array}{cc}\frac{1}{2}(\alpha I_m-P_1{\mathrm{AP}}_1^T)\hfill & -P_1{\mathrm{BP}}_2^T\hfill \\ \multicolumn{1}{c}{O}& \frac{1}{\tau }Q\hfill \end{array}\right)\left(\begin{array}{c}u\hfill \\ \multicolumn{1}{c}{v}\end{array}\right)\hfill \\ =\lambda \left(\begin{array}{cc}\frac{1}{2}(\alpha I_m+P_1{\mathrm{AP}}_1^T)\hfill & O\hfill \\ \multicolumn{1}{c}{-P_2{B}^T{P}_1^T}& \frac{1}{\tau }Q\hfill \end{array}\right)\left(\begin{array}{c}u\hfill \\ \multicolumn{1}{c}{v}\end{array}\right),\hfill \end{array}$$

or

$$\begin{array}{c}\multicolumn{1}{c}{\{\begin{array}{c}(\alpha I_m-P_1{\mathrm{AP}}_1^T)u-2P_1{\mathrm{BP}}_2^{T}v=\lambda (\alpha I_m+P_1{\mathrm{AP}}_1^T)u,\hfill \\ \multicolumn{1}{c}{\lambda \tau P_2{B}^T{P}_1^{T}u=(\lambda -1)\mathrm{Qv}.}\end{array}}\end{array}$$

(11)

• Lemma 3.2 Assume that $P_1\in {\mathbb{R}}^{m\times m}$ and $P_2\in {\mathbb{R}}^{n\times n}$ are positive definite, $A\in {\mathbb{R}}^{m\times m}$ and $Q\in {\mathbb{R}}^{n\times n}$ are symmetric and positive definite, $B\in {\mathbb{R}}^{m\times n}$ is of full column rank, $\alpha >0$ and $\tau >0$ . If $\lambda$ is an eigenvalue of $T(\alpha ,\tau )$ and $(u;v)$ is the corresponding eigenvector, then $\lambda \ne 1$ and $u\ne 0$ .

Proof. If $\lambda =1$ , then (11) becomes

$$\left(\begin{array}{cc}{P}_1{\mathrm{AP}}_1^T\hfill & P_1{\mathrm{BP}}_2^T\hfill \\ \multicolumn{1}{c}{-P_2{B}^T{P}_1^T}& O\hfill \end{array}\right)\left(\begin{array}{c}u\hfill \\ \multicolumn{1}{c}{v}\end{array}\right)=0.$$

(12)

From the assumptions, the coefficient matrix in (12) is nonsingular. Thus, $u=0$ and $v=0$ , which contradicts with the fact that $(u;v)\ne 0$ since it is an eigenvector of $T(\alpha ,\tau )$ . Therefore, it must be $\lambda \ne 1$ .

Suppose that $u=0$ . Then (11) becomes

$$\{\begin{array}{c}-P_1{\mathrm{BP}}_2^{T}v=0,\hfill \\ \multicolumn{1}{c}{(\lambda -1)\mathrm{Qv}=0.}\end{array}$$

Thus $v=0$ by $\lambda \ne 1$ , and the positive definiteness of $Q$ . This is impossible since $(u;v)\ne 0$ . Therefore, our result holds.

Now, we state and prove the convergence result of the BDP-US method.

• Theorem 3.1 Suppose that $P_1\in {\mathbb{R}}^{m\times m}$ and $P_2\in {\mathbb{R}}^{n\times n}$ are positive definite, $A\in {\mathbb{R}}^{m\times m}$ and $Q\in {\mathbb{R}}^{n\times n}$ are symmetric and positive definite, $B\in {\mathbb{R}}^{m\times n}$ is of full column rank, $\alpha >0$ and $\tau >0$ . Then the BDP-US method is convergent if the parameters $\alpha$ and $\tau$ satisfy the following condition:

$$0<\tau <\frac{2\alpha }{{\lambda }_{max}(\mathcal{B})},$$

(13)

where

$$\mathcal{B}=P_1{\mathrm{BP}}_2^T{Q}^{-1}{P}_2{B}^T{P}_1^T.$$

(14)

In particular, when $\alpha =\tau$ , the BDP-US iteration method is convergent if ${\lambda }_{max}(\mathcal{B})<2$ .

Proof. Let $\lambda$ be an eigenvalue of $T(\alpha ,\tau )$ and $(u;v)$ be the corresponding eigenvector. Then $\lambda \ne 1$ and $u\ne 0$ by Lemma 3.2, and

$$v=\frac{\lambda \tau }{\lambda -1}{Q}^{-1}{P}_2{B}^T{P}_1^{T}u$$

(15)

by the second equation of (11). In the following, we establish the relationship between the eigenvalue $\lambda$ and the parameters $\alpha$ and $\tau$ in two cases.

Case 1: $P_2{B}^T{P}_1^{T}u=0$ . Then $v=0$ by (15), and the first equation of (11) becomes

$$(\alpha I_m-P_1{\mathrm{AP}}_1^T)u-\lambda (\alpha I_m+P_1{\mathrm{AP}}_1^T)u=0.$$

From $u\ne 0$ and multiplying both sides of the above equation by $\frac{u^{*}}{u^{*}u}$ , we get

$$\alpha -\eta -\lambda (\alpha +\eta )=0,$$

where $\eta =\frac{u^{*}{P}_1{\mathrm{AP}}_1^{T}u}{u^{*}u}>0$ by the positive definiteness of $A$ . Thus, $|\lambda |=\left|\frac{\alpha -\eta }{\alpha +\eta }\right|<1$ by $\alpha >0$ .

Case 2: $P_2{B}^T{P}_1^{T}u\ne 0$ . Then substituting (15) into the first equation of (11) yields

$$(\alpha I_m+P_1{\mathrm{AP}}_1^T)u{\lambda }^2-2(\alpha I_m-\mathcal{B})u\lambda +(\alpha I_m-P_1{\mathrm{AP}}_1^T)u=0,$$

(16)

where $\mathcal{B}$ is defined in (14). Multiplying both sides of the above equation by $\frac{u^{*}}{u^{*}u}$ yields

$$(\alpha +\eta ){\lambda }^2-2(\alpha -\xi \tau )\lambda +\alpha -\eta =0,$$

(17)

where $\xi =\frac{u^{*}\mathcal{B}u}{u^{*}u}>0$ due to the positive definiteness of $Q$ .

From (17) and Lemma 3.1, we know that $|\lambda |<1$ if and only if

$$2\left|\frac{\alpha -\xi \tau }{\alpha +\eta }\right|<1+\frac{\alpha -\eta }{\alpha +\eta },$$

That is,

$$0<\tau <\frac{2\alpha }{\xi }.$$

(18)

Since $\xi \le {\lambda }_{max}(\mathcal{B})$ , (13) implies (18). Thus our results are obtained.

Spectral properties of the preconditioned matrix $\mathit{M}(\mathbf{\alpha }\mathbf{,}\mathbf{\tau }{\mathbf{)}}^{\mathbf{-}1}\hat{\mathcal{A}}$

It is worth pointing out that the BDP-US iteration method, a stationary iteration scheme, provides a preconditioner $M(\alpha ,\tau )$ defined in (8), which can accelerate the convergence rate of GMRES method when applied to the linear system (2). In this section, we shall analyze spectral properties of the preconditioned matrix $M(\alpha ,\tau {)}^{-1}\hat{\mathcal{A}}$ .

As is known, when the preconditioner $M(\alpha ,\tau )$ is applied to accelerate the convergence of Krylov subspace methods, one must solve an intermediate linear system

$$\left(\begin{array}{cc}\frac{1}{2}(\alpha I_m+P_1{\mathrm{AP}}_1^T)\hfill & O\hfill \\ \multicolumn{1}{c}{-P_2{B}^T{P}_1^T}& \frac{1}{\tau }Q\hfill \end{array}\right)\left(\begin{array}{c}{w}_1\hfill \\ \multicolumn{1}{c}{w_2}\end{array}\right)=\left(\begin{array}{c}{r}_1\hfill \\ \multicolumn{1}{c}{r_2}\end{array}\right)$$

(19)

at each iteration, where $w=(w_1^T,w_2^T{)}^T$ and $r=(r_1^T,r_2^T{)}^T\in {\mathbb{R}}^{m+n}$ . This system consists of two smaller systems and is solved as follows: First, solve the linear system

$$(\alpha I_m+P_1{\mathrm{AP}}_1^T)w_1=2r_1,$$

(20)

then solve $Q{w}_2=\tau (r_2+P_2{B}^T{P}_1^T{w}_1)$ . Noticing that $n<<m$ and $Q\in {\mathbb{R}}^{n\times n}$ , (20) can be solved by either the (sparse) Cholesky factorization or the preconditioned conjugate gradient method, and the much smaller $Q{w}_2=\tau (r_2+P_2{B}^T{P}_1^T{w}_1)$ can be solved by the Cholesky factorization.

Theorem below gives spectral properties of the preconditioned matrix $M(\alpha ,\tau {)}^{-1}\hat{\mathcal{A}}$ .

• Theorem 4.1 Assume that $P_1\in {\mathbb{R}}^{m\times m}$ and $P_2\in {\mathbb{R}}^{n\times n}$ are positive definite, $A\in {\mathbb{R}}^{m\times m}$ and $Q\in {\mathbb{R}}^{n\times n}$ are symmetric and positive definite, $B\in {\mathbb{R}}^{m\times n}$ is of full column rank. Then the real part of any eigenvalue $\hat{\lambda }$ of $M(\alpha ,\tau {)}^{-1}\hat{\mathcal{A}}$ is positive, that is, the preconditioned matrix $M(\alpha ,\tau {)}^{-1}\hat{\mathcal{A}}$ is positive stable. Furthermore, $\rho (M(\alpha ,\tau {)}^{-1}\hat{\mathcal{A}})<1$ if $\alpha >max\{{\lambda }_{max}(P_1{\mathrm{AP}}_1^T),2\tau {\lambda }_{max}(\mathcal{B})-{\lambda }_{min}(P_1{\mathrm{AP}}_1^T)\}$ , where $\mathcal{B}$ is defined by (14).

Proof. Let $\hat{\lambda }$ be an eigenvalue of $M(\alpha ,\tau {)}^{-1}\hat{\mathcal{A}}$ and $z=(u;v)$ be the corresponding eigenvector. Then

$$T(\alpha ,\tau )z=(1-\hat{\lambda })z$$

by (9) and (10). This implies that $1-\hat{\lambda }$ is an eigenvalue of $T(\alpha ,\tau )$ , and (17) holds at $\lambda =1-\hat{\lambda }$ by the proof of Theorem 3.1, that is,

$$(\alpha +\eta ){\hat{\lambda }}^2-2(\eta +\tau \xi )\hat{\lambda }+2\tau \xi =0,$$

(21)

where $\eta =u^{*}{P}_1{\mathrm{AP}}_1^{T}u/(u^{*}u)>0$ , and $\xi =u^{*}\mathcal{B}u/(u^{*}u)\ge 0$ since $u\ne 0$ by Lemma 3.2. Thus

$$\hat{\lambda }=\frac{\eta +\tau \xi \pm \sqrt{{\eta }^2-2\alpha \tau \xi +{\tau }^2{\xi }^2}}{\alpha +\eta }.$$

(22)

This, $\alpha >0$ and $\tau >0$ imply that the real part of $\hat{\lambda }$ , $\mathrm{R}e(\hat{\lambda })=(\eta +\tau \xi )/(\alpha +\eta )>0$ when $\mathrm{\Delta }={\eta }^2-2\alpha \tau \xi +{\tau }^2{\xi }^2<0$ , and $\hat{\lambda }>0$ when $\mathrm{\Delta }\ge 0$ since $(\eta +\tau \xi {)}^2-\mathrm{\Delta }=2\tau \xi (\eta +\alpha )>0$ .

Furthermore, from (21) and Lemma 3.1, $|\hat{\lambda }|<1$ if and only if $2\tau \xi <\alpha +\eta$ and $\eta <\alpha$ , that is, $\alpha >max\{\eta ,2\tau \xi -\eta \}$ . Thus the assertion holds by ${\lambda }_{min}(P_1{\mathrm{AP}}_1^T)\le \eta \le {\lambda }_{max}(P_1{\mathrm{AP}}_1^T)$ and $\xi \le {\lambda }_{max}(\mathcal{B})$ .

The parameter choice

In general, the effectiveness of an iteration method depends largely on the choice of its parameters, and it is very important to find a good strategy for approaching the optimal parameters of the iteration method since it is very hard to find the optimal parameters. In particular, Bai ²⁸ provided the quasi-optimal parameters of HSS iteration methods to solve the large sparse saddle-point problems. Following the similar analysis as Bai, ²⁸ we shall give a strategy to select the iteration parameters of the BDP-US method for solving the linear system (2).

Let $\lambda$ be an eigenvalue of $T(\alpha ,\tau )$ and $(u;v)$ be a corresponding eigenvector. When $B$ is of full column rank, $u\ne 0$ and $\lambda \ne 1$ by Lemma 3.2, and

$$\lambda =\frac{\alpha -\tau \xi \pm \sqrt{{\eta }^2-2\alpha \tau \xi +{\tau }^2{\xi }^2}}{\alpha +\eta }$$

by (17), where $\eta =u^{*}{P}_1{\mathrm{AP}}_1^{T}u/(u^{*}u)>0$ and $\xi =u^{*}\mathcal{B}u/(u^{*}u)>0$ . In the following, we analyze the choices for the parameters in two cases.

Case 1: $\mathrm{\Delta }={\eta }^2-2\alpha \tau \xi +{\tau }^2{\xi }^2\ge 0$ . Then $\alpha <\eta$ , or $\alpha \ge \eta$ and $\tau \ge (\alpha +\sqrt{{\alpha }^2-{\eta }^2})/\xi$ , or $\alpha \ge \eta$ and $\tau \le (\alpha -\sqrt{{\alpha }^2-{\eta }^2})/\xi$ . So

$$\begin{array}{c}|\lambda |=\hfill \\ \{\begin{array}{cc}\frac{\alpha -\tau \xi +\sqrt{\mathrm{\Delta }}}{\alpha +\eta },\hfill & \mathrm{if}\alpha \ge max\{\eta ,\tau \xi +\sqrt{{\alpha }^2-{\eta }^2}\},\mathrm{or}\eta >\alpha \ge \tau \xi ,\hfill \\ \multicolumn{1}{c}{\frac{\tau \xi -\alpha +\sqrt{\mathrm{\Delta }}}{\alpha +\eta },}& \mathrm{if}\eta \le \alpha \le \tau \xi -\sqrt{{\alpha }^2-{\eta }^2},\mathrm{or}\alpha \le min\{\eta ,\tau \xi \}.\hfill \end{array}\hfill \end{array}$$

• a)  $\alpha \ge max\{\eta ,\tau \xi +\sqrt{{\alpha }^2-{\eta }^2}\}$ . Then $|\lambda |\le (\alpha -\tau \xi )/\eta$ by $\alpha \ge \eta$ and $\mathrm{\Delta }\le (\alpha -\tau \xi {)}^2$ . Due to the linearity of $\alpha -\tau \xi$ with respect to $\alpha$ and $\tau$ , $minmax\{|\lambda |\}=0$ at $\alpha =\tau \xi =\eta$ when $\alpha \ge \eta \ge \tau \xi +\sqrt{{\alpha }^2-{\eta }^2}$ , and $minmax\{|\lambda |\}=\sqrt{{(\alpha /\eta )}^2-1}$ at $\alpha =\tau \xi +\sqrt{{\alpha }^2-{\eta }^2}$ when $\alpha \ge \tau \xi +\sqrt{{\alpha }^2-{\eta }^2}\ge \eta$ .

• b)  $\eta >\alpha \ge \tau \xi$ . Then

$$|\lambda |\le g(\alpha ,\tau )=\sqrt{\frac{\eta -\tau \xi }{\eta +\tau \xi }}+\frac{\alpha -\tau \xi }{\alpha +\tau \xi }$$

by $\mathrm{\Delta }\le {\eta }^2-(\tau \xi {)}^2$ since $-\alpha \tau \xi \le -(\tau \xi {)}^2$ . Since

$$\{\begin{array}{c}\frac{\partial g(\alpha ,\tau )}{\partial \alpha }=\frac{2\tau \xi }{{(\alpha +\tau \xi )}^2}>0,\hfill \\ \multicolumn{1}{c}{\frac{\partial g(\alpha ,\tau )}{\partial \tau }=-\frac{2\alpha \xi }{{(\alpha +\tau \xi )}^2}-\frac{\xi \eta }{\sqrt{{(\eta +\tau \xi )}^3(\eta -\tau \xi )}}<0,}\end{array}$$

$g(\alpha ,\tau )$ is increasing and decreasing with respect to $\alpha$ and $\tau$ , respectively. Thus, $g(·,\tau )$ and $g(\alpha ,·)$ attain the minimum when $\alpha =\tau \xi$ , which leads to $minmax\{|\lambda |\}=\sqrt{(\eta -\alpha )/(\eta +\alpha )}$ .

• c)  $\eta \le \alpha \le \tau \xi -\sqrt{{\alpha }^2-{\eta }^2}$ . Then $|\lambda |\le (\tau \xi -\alpha )/\eta$ , and $minmax\{|\lambda |\}=\sqrt{{(\alpha /\eta )}^2-1}$ at $\alpha =\tau \xi -\sqrt{{\alpha }^2-{\eta }^2}$ .

• d)  $\alpha \le min\{\eta ,\tau \xi \}$ . Let $f(\alpha ,\tau )=(\tau \xi -\alpha +\sqrt{\mathrm{\Delta }})/(\alpha +\eta )$ . Then

$$\{\begin{array}{c}\frac{\partial f(\alpha ,\tau )}{\partial \alpha }=-\frac{\xi \tau +\sqrt{\mathrm{\Delta }}}{(\alpha +\eta )\sqrt{\mathrm{\Delta }}}-\frac{\xi \tau -\alpha +\sqrt{\mathrm{\Delta }}}{{(\alpha +\eta )}^2}<0,\hfill \\ \multicolumn{1}{c}{\frac{\partial f(\alpha ,\tau )}{\partial \tau }=\frac{\xi \tau -\alpha +\sqrt{\mathrm{\Delta }}}{(\alpha +\eta )\sqrt{\mathrm{\Delta }}}\xi >0.}\end{array}$$

This implies that $f(\alpha ,\tau )$ is decreasing and increasing with respect to $\alpha$ and $\tau$ , respectively. Thus, $minmax\{|\lambda |\}=0$ at $\alpha =\tau \xi =\eta$ when $\tau \xi \ge \eta \ge \alpha$ ; and $minmax\{|\lambda |\}=\sqrt{(\eta -\alpha )/(\eta +\alpha )}$ at $\alpha =\tau \xi$ when $\eta \ge \tau \xi \ge \alpha$ .

Case 2: $\mathrm{\Delta }={\eta }^2-2\alpha \tau \xi +{\tau }^2{\xi }^2<0$ . That is, ${\eta }^2+{\tau }^2{\xi }^2<2\alpha \tau \xi$ . Then

$$|\lambda |=\left|\sqrt{{\left(\frac{\alpha -\tau \xi }{\alpha +\eta }\right)}^2+\frac{2\alpha \tau \xi -{\eta }^2-{\tau }^2{\xi }^2}{{(\alpha +\eta )}^2}}\right|=\left|\sqrt{\frac{\alpha -\eta }{\alpha +\eta }}\right|<1.$$

In summary, $minmax\{|\lambda |\}=\sqrt{{(\alpha /\eta )}^2-1}$ at $\alpha =\tau \xi \pm \sqrt{{\alpha }^2-{\eta }^2}\ge \eta$ with $\alpha <\sqrt{2}\eta$ , and $minmax\{|\lambda |\}=\sqrt{(\eta -\alpha )/(\eta +\alpha )}$ at $\alpha =\tau \xi$ . Therefore, the best choices of $\alpha$ and $\tau$ are $\alpha =\eta =\tau \xi$ .

Numerical experiments

In this section, a numerical example is provided to support the theoretical results and illustrate the effectiveness of BDP-US method. All experiments are conducted in MATLAB (R2015a) and terminated if the current residual

$$\mathrm{RES}:=\frac{\sqrt{\parallel P_{1}f-P_{1}A{x}^k-P_{1}B{y}^k{\parallel }_2^2+\parallel P_{2}g-P_2{B}^T{x}^k{\parallel }_2^2}}{\sqrt{\parallel P_{1}f{\parallel }_2^2+\parallel P_{2}g{\parallel }_2^2}}$$

is less than ${10}^{-7}$ or the number of iteration steps is greater than the prescribed number $k_{max}=1000$ . In our computations, the initial guess is set to be zero, and the right-hand side $\hat{b}$ is chosen such that the exact solution of the saddle-point problem (2) is $(1,1,1,...,1{)}^T\in {\mathbb{R}}^{m+n}$ .

• Example 6.1 ¹⁸ Consider the following linear system, where

$$\begin{array}{c}A=\left(\begin{array}{cc}{I}_p\otimes T+T\otimes I_p\hfill & O\hfill \\ \multicolumn{1}{c}{O}& I_p\otimes T+T\otimes I_p\hfill \end{array}\right)\in {\mathbb{R}}^{2p^2\times 2p^2}\mathrm{and}\hfill \\ B=\left(\begin{array}{c}{I}_p\otimes F\hfill \\ \multicolumn{1}{c}{F\otimes I_p}\end{array}\right)\in {\mathbb{R}}^{2p^2\times p^2},\hfill \end{array}$$

$T=\mathrm{tridiag}(-1,2,-1)/h^2\in {\mathbb{R}}^{p\times p},F=\mathrm{tridiag}(-1,1,0)/h\in {\mathbb{R}}^{p\times p},$ $h=1/(p+1)$ is the discretization mesh-size, ⊗ is the Kronecker product symbol, $m=2p^2$ and $n=p^2$ .

To show the advantages of BDP-US method and preconditioner $M(\alpha ,\tau )$ , we compare it with three iteration methods (Uzawa-SOR, ²⁶ Uzawa-Low, ²⁹ and one-parameter variant of preconditioned Uzawa (OVPU) ³⁰ methods) and their corresponding preconditioners. In all experiments, the parameter $\omega$ in OVPU method is optimal value, while the parameters in the other methods are set to the experimentally computed optimal ones, which leads to the least number of iteration steps.

Take $P_1$ , $P_2$ , and $Q$ as shown in Table 1, $p=8,16$ and $24$ . Tables 2–4 report the numerical results obtained by these methods for two cases, where “IT” denotes the number of iteration steps and “CPU” is the elapsed CPU time in seconds. From Tables 2 to 4, we see that BDP-US method outperforms the other three methods, and BDP-US-I has a worse performance than BDP-US-II. Thus, the BDP-US method with suitable parameters $\alpha ,\tau$ , and $Q$ is more efficient.

Table 1.

Choices of the matrices of $P_1$ , $P_2$ , and $Q$ .

Cases	$P_1$	$P_2$	$Q$
I	$I_m$	$I_n$	$B^T\mathrm{tridiag}(A{)}^{-1}B$
II	$\mathrm{tridiag}(A{)}^{-1}$	$I_n$	$B^T\mathrm{tridiag}(A{)}^{-1}B$

Table 2.

Numerical results for example 6.1 ( $p=8$ ).

Cases	$(\omega ,\alpha ,\tau )$	IT	CPU	RES
Uzawa-SOR-I	(1.38,–,0.10)	159	0.0313	9.87E-08
Uzawa-SOR-II	(1.20,–,0.17)	128	0.0288	9.85E-08
Uzawa-Low-I	(1.50,–,0.13)	139	0.0307	8.88E-08
Uzawa-Low-II	(1.52,–,0.26)	87	0.0199	9.90E-08
OVPU-I	(0.4136,–,–)	81	0.0192	9.60E-08
OVPU-II	(0.4136,–,–)	79	0.0190	9.99E-08
BDP-US-I	(–,118,0.248)	68	0.0187	9.83E-08
BDP-US-II	(–,0.50,0.50)	49	0.0084	8.39E-08

Table 3.

Numerical results for example 6.1 ( $p=16$ ).

Cases	$(\omega ,\alpha ,\tau )$	IT	CPU	RES
Uzawa-SOR-I	(1.41,–,0.032)	436	3.5673	9.80E-08
Uzawa-SOR-II	(1.25,–,0.056)	366	2.1211	9.69E-08
Uzawa-Low-I	(1.48,–,0.038)	418	2.3622	9.73E-08
Uzawa-Low-II	(1.48,–,0.075)	272	2.3017	9.78E-08
OVPU-I	(0.1529,–,–)	277	1.7060	9.82E-08
OVPU-II	(0.1529,–,–)	266	1.6631	9.27E-08
BDP-US-I	(–,237,0.137)	112	0.6171	9.30E-08
BDP-US-II	(–,0.308,0.308)	86	0.4954	7.90E-08

Table 4.

Numerical results for example 6.1 ( $p=24$ ).

Cases	$(\omega ,\alpha ,\tau )$	IT	CPU	RES
Uzawa-SOR–I	(1.45,–,0.015)	831	22.8820	9.89E-08
Uzawa-SOR–II	(1.23,–,0.028)	711	18.5213	9.92E-08
Uzawa-Low–I	(1.46,–,0.019)	809	22.6989	9.22E-08
Uzawa-Low–II	(1.45,–,0.041)	520	17.2806	9.11E-08
OVPU-I	(0.0764,–,–)	602	14.0528	9.82E-08
OVPU-II	(0.0764,–,–)	569	12.4945	9.73E-08
BDP-US-I	(–,348,0.093)	154	5.1158	9.48E-08
BDP-US-II	(–,0.238,0.238)	114	2.7059	8.95E-08

Moreover, we accelerate the convergence rate of the GMRES(10) method by using the $M(\alpha ,\tau )$ , Uzawa-SOR, Uzawa-Low, and OVPU preconditioners (denote by ${\mathcal{P}}_1$ , ${\mathcal{P}}_2$ , ${\mathcal{P}}_3$ , and ${\mathcal{P}}_4$ respectively). Tables 5 and 6 list the outer(inner) iteration numbers IT, the CPU and the RES for accelerating these iteration methods by GMRES(10) in two cases when $p=16$ , where IT, CPU, and RES are the same as Table 2. From Tables 5 and 6, we see that GMRES(10) with preconditioner ${\mathcal{P}}_1$ requires less IT and CPU than that with the other three compared preconditioners, and Case II outperforms Case I. Thus, GMRES with the preconditioner ${\mathcal{P}}_1$ achieves faster convergence rate.

Table 5.

Outer(inner) iterations, CPU(s), and RES of GMRES(10) for Case I when $p=16$ .

Preconditioners	IT	CPU	RES
Uzawa-SOR $({\mathcal{P}}_2)$	13(9)	2.7271	9.46E-08
Uzawa-Low $({\mathcal{P}}_3)$	13(7)	2.5116	9.34E-08
OVPU $({\mathcal{P}}_4)$	18(9)	3.2805	7.94E-08
BDP-US $({\mathcal{P}}_1)$	5(3)	0.8190	7.93E-08

Table 6.

Outer(inner) iterations, CPU(s), and RES of GMRES(10) for Case II when $p=16$ .

Preconditioners	IT	CPU	RES
Uzawa-SOR $({\mathcal{P}}_2)$	12(8)	2.3917	9.50E-08
Uzawa-Low $({\mathcal{P}}_3)$	11(9)	2.2694	9.36E-08
OVPU $({\mathcal{P}}_4)$	21(8)	3.8201	8.01E-08
BDP-US $({\mathcal{P}}_1)$	4(7)	0.6848	7.72E-08

To see more clearly, Figures 1 and 2 depict the eigenvalue distributions of the corresponding preconditioned matrices for two cases when $p=16$ . Clearly, the real parts of $\lambda ({\mathcal{P}}_1^{-1}\hat{\mathcal{A}})$ are all nonnegative, which coincides with that of Theorem 4.1, and the eigenvalues of ${\mathcal{P}}_1^{-1}\hat{\mathcal{A}}$ are more clustered than the other three preconditioned matrices. Therefore, the preconditioner ${\mathcal{P}}_1$ with suitable parameters is more effective than ${\mathcal{P}}_2$ , ${\mathcal{P}}_3$ , and ${\mathcal{P}}_4$ when applied to solve large sparse saddle-point problems.

Figure 1.

Eigenvalue distributions of the corresponding preconditioned matrices for Case I when $p=16$ : (a) Uzawa-SOR, (b) Uzawa-Low, (c) OVPU, and (d) BDP-US.

Figure 2.

Eigenvalue distributions of the corresponding preconditioned matrices for Case II when $p=16$ : (a) Uzawa-SOR, (b) Uzawa-Low, (c) OVPU, and (d) BDP-US.

Conclusions

This paper develops a block-diagonally preconditioned US method for solving the saddle-point problems by generalizing the US method in Li and Ma ²⁵ and sufficiently utilizing the block diagonal preconditioning technique, and provides a sufficient condition to ensure its convergence. Then, a splitting preconditioner is presented, the spectral properties of the preconditioned matrix are analyzed, and the choice of the parameters for this iteration method is discussed. Finally, numerical experiments are conducted to demonstrate that BDP-US method performs better than some existing methods and the preconditioner $M(\alpha ,\tau )$ can accelerate the corresponding Krylov subspace methods. Moreover, further explorations should be focused on finding optimal parameters of the BDP-US method and giving more accurate description for the spectral distribution of the preconditioned matrix.