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

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.


Introduction
Consider the linear system of the 2 3 2 block structure: where A 2 R m 3 m is symmetric and positive definite, B 2 R m 3 n with n ( m, x, f 2 R m , y, g 2 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., 1 Bai and Bai, 2 Bai and Pan 3 ).
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][5][6][7][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][10][11][12] SOR (successive overrelaxation) type methods, [13][14][15][16] HSS (Hermitian and skew-Hermitian splitting) type methods and its accelerated variants, [17][18][19][20]28 and Krylov subspace methods with suitable preconditioners. [21][22][23][24] In particular, Bai 28 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 21 and the classical Uzawa method, 9 Li and Ma 25 proposed the following Uzawa splitting (US) iteration method x k + 1 = x k + 2(aI m + A) À1 ( f À Ax k À By k ) y k + 1 = y k + tQ À1 (B T x k + 1 À g) to solve singular saddle-point problem (1) with rank(B) ł n and rank(Á) being the rank of the matrix, where Q 2 R n 3 n is a symmetric positive definite matrix, a and t are two positive constants. Even though this method is a special case of the parameterized inexact Uzawa (PIU) method, 15 it has better performance than the Uzawa method and Uzawa-SOR methods. 26 Moreover, this method is convergent when where l max is the largest eigenvalue of the matrix BQ À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 l max may become larger, which makes a too large, or t too small. For example, when p = 8, 16, 24, the given values of (a, t) in Li and Ma 25 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, 25 and it requires to accurately solve a linear system with the coefficient matrix aI 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 25 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 Re(Á) the real part of a complex number, by the superscripts T the transpose of a matrix or vector, by r(Á) the spectral radius of a matrix, by (u; v) the vector (u Ã , v Ã ) Ã and by l 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 where P 1 2 R m 3 m and P 2 2 R n 3 n are positive definite matrices. Then by left and right multiplying the linear system (1) with P À1 and P ÀT respectively, we obtain which can be written aŝ ð2Þ Similar to Bai et al., 21 we can define the corresponding shift-splitting of P 1 AP T 1 as where a is a positive constant. For the available approximationẑ k = ((P ÀT 1 x k ) T , (P ÀT 2 y k ) T ) T , k = 0, 1, . . ., the following iteration scheme is first applied to the first m equations in (2) to compute the update x k + 1 : Then the Uzawa iteration method is used to the last n equations in (2) to compute the update y k + 1 : where Q 2 R n 3 n is symmetric and positive definite, and t is a positive constant. For simplicity, we denotex = P ÀT 1 x,ŷ = P ÀT 2 y, x k = P ÀT 1 x k andŷ k = P ÀT 2 y k . Then the results in (4) and (5) suggest the following BDP-US method for the linear system (2): where the iteration matrix On the other hand, let and then we obtain a splitting A = M(a, t) À N (a, t), and Overall, we can perform the BDP-US method as follows: Given an initial guessẑ 0 = ((x 0 ) T , (ŷ 0 ) T ) T witĥ x 0 2 R m andŷ 0 2 R n , two positive parameter a and t; for k = 0, 1, 2, . . ., computeẑ k + 1 by the iteration scheme (6) until the iteration sequence fẑ k g ' converges. Thus we obtain the converged sequence Clearly, this BDP-US method reduces to the US method in Li and Ma 25 by setting P 1 = I m and P 2 = I n , and different Q, a and t 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. It is well known that the BDP-US method (7) is convergent if and only if r(T(a, t))\1. Let l be an eigenvalue of T (a, t), and (u; v) with u 2 C m and v 2 C n be the corresponding eigenvector, then This, (8) and (9) imply that Lemma 3.2 Assume that P 1 2 R m 3 m and P 2 2 R n 3 n are positive definite, A 2 R m 3 m and Q 2 R n 3 n are symmetric and positive definite, B 2 R m 3 n is of full column rank, a.0 and t.0. If l is an eigenvalue of T(a, t) and (u; v) is the corresponding eigenvector, then l 6 ¼ 1 and u 6 ¼ 0.
Proof. If l = 1, then (11) becomes From the assumptions, the coefficient matrix in (12) is nonsingular. Thus, u = 0 and v = 0, which contradicts with the fact that (u; v) 6 ¼ 0 since it is an eigenvector of T (a, t). Therefore, it must be l 6 ¼ 1.
Suppose that u = 0. Then (11) becomes Thus v = 0 by l 6 ¼ 1, and the positive definiteness of Q. This is impossible since (u; v) 6 ¼ 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 2 R m 3 m and P 2 2 R n 3 n are positive definite, A 2 R m 3 m and Q 2 R n 3 n are symmetric and positive definite, B 2 R m 3 n is of full column rank, a.0 and t.0. Then the BDP-US method is convergent if the parameters a and t satisfy the following condition: where In particular, when a = t, the BDP-US iteration method is convergent if l max (B)\2.
Proof. Let l be an eigenvalue of T (a, t) and (u; v) be the corresponding eigenvector. Then l 6 ¼ 1 and u 6 ¼ 0 by by the second equation of (11). In the following, we establish the relationship between the eigenvalue l and the parameters a and t in two cases. Case 1: (15), and the first equation of (11) becomes (aI m À P 1 AP T 1 )u À l(aI m + P 1 AP T 1 )u = 0: From u 6 ¼ 0 and multiplying both sides of the above equation by u Ã u Ã u , we get .0 by the positive definiteness of A. Thus, jlj = a À h a + h \1 by a.0.
Case 2: P 2 B T P T 1 u 6 ¼ 0. Then substituting (15) into the first equation of (11) yields where B is defined in (14). Multiplying both sides of the above equation by where j = u Ã Bu u Ã u .0 due to the positive definiteness of Q.
From (17) and Lemma 3.1, we know that jlj\1 if and only if Since j ł l max (B), (13) implies (18). Thus our results are obtained.
Spectral properties of the preconditioned matrix M(a, t) À1Â It is worth pointing out that the BDP-US iteration method, a stationary iteration scheme, provides a preconditioner M(a, t) 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(a, t) À1Â . As is known, when the preconditioner M(a, t) is applied to accelerate the convergence of Krylov subspace methods, one must solve an intermediate linear system at each iteration, where w = (w T 1 , w T 2 ) T and r = (r T 1 , r T 2 ) T 2 R m + n . This system consists of two smaller systems and is solved as follows: First, solve the linear system then solve Qw 2 = t(r 2 + P 2 B T P T 1 w 1 ). Noticing that n ( m and Q 2 R n 3 n , (20) can be solved by either the (sparse) Cholesky factorization or the preconditioned conjugate gradient method, and the much smaller Qw 2 = t(r 2 + P 2 B T P T 1 w 1 ) can be solved by the Cholesky factorization.
Theorem below gives spectral properties of the preconditioned matrix M(a, t) À1Â .
Theorem 4.1 Assume that P 1 2 R m 3 m and P 2 2 R n 3 n are positive definite, A 2 R m 3 m and Q 2 R n 3 n are symmetric and positive definite, B 2 R m 3 n is of full column rank. Then the real part of any eigenvaluel of M(a, t) À1Â is positive, that is, the preconditioned matrix M(a, t) À1Â is positive stable. Furthermore, r(M(a, t) À1Â )\1 if a. maxfl max (P 1 AP T 1 ), 2tl max (B) À l min (P 1 AP T 1 )g, where B is defined by (14).

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 28 provided the quasi-optimal parameters of HSS iteration methods to solve the large sparse saddlepoint problems. Following the similar analysis as Bai, 28 we shall give a strategy to select the iteration parameters of the BDP-US method for solving the linear system (2). Let l be an eigenvalue of T (a, t) and (u; v) be a corresponding eigenvector. When B is of full column rank, u 6 ¼ 0 and l 6 ¼ 1 by Lemma 3.2, and l = a À tj 6 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi h 2 À 2atj + t 2 j 2 p a + h by (17), where h = u Ã P 1 AP T 1 u=(u Ã u).0 and j = u Ã Bu= (u Ã u).0. In the following, we analyze the choices for the parameters in two cases.
g(a, t) is increasing and decreasing with respect to a and t, respectively. Thus, g( Á , t) and g(a, Á ) attain the minimum when a = tj, which leads to . Then jlj ł (tj À a)=h, Then ∂f (a, t) This implies that f (a, t) is decreasing and increasing with respect to a and t, respectively. Thus, min maxfjljg = 0 at a = tj = h when tj ø h ø a; and min maxfjljg = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi (h À a)=(h + a) p at a = tj when h ø tj ø a.

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 p 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 sideb is chosen such that the exact solution of the saddle-point problem (2) is (1, 1, 1, :::, 1) T 2 R m + n .
Example 6.1 18 Consider the following linear system, where 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(a, t), we compare it with three iteration methods (Uzawa-SOR, 26 Uzawa-Low, 29 and oneparameter variant of preconditioned Uzawa (OVPU) 30 methods) and their corresponding preconditioners. In all experiments, the parameter v 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 a, t, and Q is more efficient.
Moreover, we accelerate the convergence rate of the GMRES(10) method by using the M(a, t), Uzawa-SOR, Uzawa-Low, and OVPU preconditioners (denote by P 1 , P 2 , P 3 , and 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 P 1 requires less IT and CPU than that with the other three  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 l(P À1 1Â ) are all nonnegative, which coincides with that of Theorem 4.1, and the eigenvalues of P À1 1Â are more clustered than the other three preconditioned matrices. Therefore, the preconditioner P 1 with suitable   parameters is more effective than P 2 , P 3 , and P 4 when applied to solve large sparse saddle-point problems.

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 25 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(a, t) 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.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.