We consider the saddle point problems with the following block three-by-three structure
where is a symmetric positive definite matrix, and are matrices of full row rank, , and are given vectors. The linear system of equations (1) has many practical applications in scientific computing and engineering fields, including the constrained quadratic optimization problem, the constrained least squares problem and the computational fluid dynamics .
In practice, the coefficient matrix in (1) is generally sparse and of large size. Therefore, iterative methods are usually preferred for their low memory requirement per iteration. Among these candidates, the Uzawa-type methods , the matrix splitting methods  and the Krylov subspace methods  have received much attention during the past decades. For instance, the Uzawa method  is a classical approach which is known for its low computational overhead and easy implementation. To improve the numerical performance, some Uzawa-type variants have been proposed recently which include the inexact Uzawa   , Uzawa-SOR  and Uzawa-PSS , etc. The other way for solving the saddle point problems is to exploit the matrix splitting. By using the Hermitian and skew-Hermitian splitting (HSS), Bai et al.  propose the HSS iteration method. Following this reasoning, some efficient variants such as the accelerated Hermitian and skew-Hermitian splitting (AHSS) , the modified Hermitian and skew-Hermitian splitting (MHSS)  and the relaxed Hermitian and skew- Hermitian splitting (RHSS)  have been proposed henceforth. Recently, the sophisticated Krylov subspace methods have been used to solve the saddle point problems. Due to the structure of the saddle point matrices, however, a slow convergence or even stagnation is often witnessed during the actual computation. Such drawback can be alleviated by employing an appropriate preconditioner. In the context of saddle point problems, many efficient preconditioners, to name just a few, the block diagonal preconditioners, the constrained preconditioners, the HSS-based preconditioners, have been developed; see     and the more recent monograph  for detail.
In , Bai et al. propose a class of shift-splitting iteration method for solving non-Hermitian positive definite linear system . It can be recast in the following way. Suppose that is a non-Hermitian positive definite matrix. Then the shift-splitting of the coefficient matrix A is given by
which leads naturally to the shift-splitting iteration scheme
Theoretical analysis shows that shift-splitting iteration is endowed with an unconditional convergence. By using this sound convergence property, Cao et al. come up with the shift-splitting to solve the saddle point problems and propose a shift-splitting preconditioner and a local shift-splitting variant . Numerical experiments show that the proposed two shift-splitting-type preconditioners are superior to the HSS preconditioner in terms of CPU time and iteration steps. The shift-splitting technique has also been generalized by incorporating two relaxation parameters . Furthermore, Salkuyeh et al. extend it to solve the generalized saddle point problems and propose a modified generalized shift-splitting preconditioner .
By taking into account the special structure of (1), Cao  proposes shift- splitting preconditioners for solving (1). Motivated by  and , we devise a new generalized shift-splitting (GSS) iteration scheme for the structured linear system (1). A corresponding GSS preconditioner is induced to accelerate the convergence of GMRES. Theoretical results on the convergence of GSS iteration method and the eigenvalue distribution of the preconditioned matrix are also given.
The organization of this paper is as follows. In Section 2, we introduce the generalized shift-splitting iteration method for solving the block three-by-three saddle point problem and look into its convergence property. In Section 3, we obtain a GSS preconditioner derived from the GSS iteration and analyze the spectrum of the preconditioned matrix. In Section 4, some numerical experiments are carried out to illustrate the effectiveness of GSS by comparing it with other shift-splitting preconditioners given in . Finally, some conclusions are given in Section 5.
2. Generalized Shift-Splitting Iteration Method and Its Convergence
The proposed generalized shift-splitting (GSS) iteration is based on the work of . Therefore, it is reasonable for us to recap briefly the shift splitting used in  before introducing the new iteration method in this section.
In its simplest form, the shift splitting (SS) of the coefficient matrix in (1) goes as follows
where I is the identity matrix with conforming size. The parameter in (2) is used to monitor the convergence of the SS iteration method and hence the spectrum of the preconditioned matrix. In other words, plays the same role in adjusting the two block diagonal matrices in (2), i.e.,
However, the reason for choosing a single parameter to control simultaneously these two block diagonal matrices above remains elusive. Instead, it seems more appealing to investigate the convergence if the two block diagonal matrices are affected by two different parameters separately. This is the main motivation of our work. To this end, we refine the splitting (2) by embedding an extra parameter in the -block, which yields the following generalized shift- splitting
Here , and . Then it is trivial to find that the generalized shift-splitting (4) reduces to the shift-splitting (2) if .
Using the splitting (4), we obtain the following stationary iteration
where , and . To put it another way, the generalized shift-splitting method (GSS) for the block three-by- three saddle point problem (1) that can be stated as follows. Given an initial vector , compute
until converges, where .
It is known that the GSS iteration method (5) enjoys the unconditional convergence if the spectral radius of the iteration matrix is less than 1. The following results are adapted from  and  which are instrumental in proving such sound convergence property.
Lemma 1. Suppose that is a symmetric positive definite matrix, and are of full row rank. Let be an eigenvalue of the saddle point matrix , then
Proof. The proof can be found in ( , Lemma 2.2).
Lemma 2. Let be a non-Hermitian positive semidefinite matrix and be the extrapolated Cayley transform with . Then
1) the spectral radius of is bounded by 1, i.e., ;
2) if and only if the matrix W has a (reducing) eigenvalue of the form with and the imaginary unit.
Proof. See ( , Theorem 2.1) for detail.
The following theorem on the convergence of GSS iteration method comes out as a result of Lemma 1 and Lemma 2.
Theorem 1. Suppose that is a symmetric positive definite matrix, and are matrices of full row rank. Let and . Then the spectral radius of defined in (5) satisfies
i.e., the GSS iteration method converges to the unique solution of (1) unconditionally.
Proof. By (4) and (5), we rewrite the iteration matrix as
Let . It is easy to see that is similar to the matrix . Therefore, it suffices to prove . By the construction of , we have
Since and are positive, then is a symmetric positive definite matrix, and are matrices of full row rank. Thus is a
nonsymmetric positive semi-definite matrix. Moreover, it follows from Lemma 1 that . Therefore, it can be deduced from Lemma 2 that . The equality holds if and only if has a (reducing) eigenvalue of the form , which is impossible because . It turns out that . As a result, which completes the proof.
3. GSS Preconditioner and Its Properties
Though the GSS iteration method is proved to converge unconditionally, its convergence can be rather slow when compared with the involved Krylov subspace methods . In the context of linear system solvers, it is more interesting to exploit the generalized shift-splitting as a preconditioner to accelerate the Krylov subspace methods like GMRES.
In this section, we present the following generalized shift-splitting preconditioner (GSS)
When applying a preconditioner, the convergence of the preconditioned linear system is closely related with the eigenvalue distribution. The following finding shows that the eigenvalues of GSS-preconditioned matrix cluster in a circle centered at the point with radius less than 1. We adopt the right preconditioning here such that the residual of the original saddle point problem is the same as that of its right-preconditioned counterpart.
Theorem 2. Suppose that is a symmetric positive definite matrix, and are matrices of full row rank. Let be an eigenvalue of , then
Proof. Since the preconditioned matrix is similar to , then it is only necessary to investigate the spectrum of . Since , then it follows from (4) that
From Theorem 1, we know that , where is an eigenvalue of .
Before ending this section, we give some implementation details when applying the GSS preconditioner. In actual computations, a matrix-vector product is often needed which can be solved from the linear system . Note that the preconditioner can be decomposed as
Let and , where , and . By using the decomposition (8), we obtain the following subroutine for solving .
1) Solve from ;
2) Solve from ;
3) Solve from ;
4) Update ;
5) Update .
In computing , the computational overhead by GSS is the same as that by SS . As indicated in the numerical examples, however, an appropriate choice of generally yields better numerical performance (concerning CPU time and iteration steps) than using a single parameter (as in SS). In this sense, the GSS can be regarded as a more promising variant than its predecessor SS.
4. Numerical Experiments
In this section, we give some numerical experiments to illustrate the effectiveness of the new preconditioner GSS in terms of the number of iteration steps (iter), CPU time in seconds (time) and the relative residual norm (res). For completeness, we display both the number of restarts and the number of inner steps in the last restart; see, for instance, Table 3. The unconditioned GMRES(k), the shift-spliting-preconditioned GMRES(k) (SS)  and the relaxed shift- splitting-preconditioned GMRES(k) (RSS)  are employed for comparison, where k is the restarting frequency. In what follows, the right preconditioning is used for SS, RSS and GSS, the reason of which is given in Section 3. We set the restarting frequency k to be 5 and the initial guess to be the zero vector. The right-hand side vector b in (1) is given by . All algorithms mentioned above terminate if within 1500 restarts, where is the ith iterate. Numerical experiments are carried out on a laptop using MATLAB 2014b with Intel Core (8G RAM) under Windows 10 system.
The testing problem is adapted from . It is a block three-by-three saddle point problem with A, B and C of the following structure
where , , with being the Kronecker product and the grid size.
4.1. Choice of the Parameter Pair
As mentioned in Section 2, the main motivation for embracing an extra parameter in GSS is to strike a balance in controlling the sub-matrices in (3). A proper choice of and will make the GSS preconditioner more appealing. As shown in (7), the GSS preconditioner becomes closer to (half of) the saddle point matrix as and decrease. This implies that the choice seems to be the optimal one. However, we do not suggest using extremely small values of and for numerical concerns. Instead, we use the experimental (optimal) values of and in the following experiments.
To this end, we first tabulate the CPU time of the GSS preconditioner with varying and for and in Table 1 and Table 2, respectively. The parameter pair yields the shortest CPU time in both tables. To further verify this finding, we carry out more experiments in Figure 1 and Figure 2 by depicting the number of total iteration steps1 and CPU time of GSS with varying when is chosen to be 0.01, 0.1 and 1, respectively. Some remarks are in order. Loosely speaking, for a fixed value of , the number of total iteration steps and CPU time increase as goes up although there are some zigzags; see Figure 1 and Figure 2. Moreover, a similar pattern is also observed for when the value of is fixed while varies; for instance, when , the total iteration steps increase as changes from 0.01 to 1, which reads readily by comparing the subplots (a), (c) and (e) in Figure 1. This justifies our choice of using as the experimental optimal value.
Based on the above experiments, we use as the experimental optimal choice for and in GSS. As for the optimal value of
Table 1. CPU time of the GSS preconditioner with different values of and ( ).
in SS and RSS, we follow the choice used in  and regard it as the empirical optimal value.
4.2. Comparison with Other SS-Type Methods
In this subsection, we compare the numerical performance of GSS with SS and
Table 2. CPU time of the GSS preconditioner with different values of and ( ).
Figure 1. Total iteration steps (left) and CPU time (right) of the GSS preconditioner for different values of and ( ).
Figure 2. Total iteration steps (left) and CPU time (right) of the GSS preconditioner for different values of and ( ).
RSS regarding CPU time and iteration steps. It should be noted that all sub-linear equations involved in computing matrix-vector product are solved by direct methods.
It is known that the spectrum imposes a great effect on the convergence of the preconditioned linear system. To begin with the discussion, we plot respectively the spectral distributions of the original saddle point matrix , the GSS-pre- conditioned matrix , the SS-preconditioned matrix , and the RSS-preconditioned matrix in Figure 3. As illustrated, the eigenvalues of the original saddle point matrix are well separated, which makes the unpreconditioned GMRES method cumbersome. After employing the SS-type preconditioners, however, the eigenvalues of the preconditioned matrices become much more clustered. In particular, the eigenvalues of the GSS-precondi- tioned matrix are more clustered than those of SS and RSS; see Figure 3. Thus it can be expected that GSS will bring about an improvement over SS and RSS in solving (1). This intuition is echoed by Table 3, where the saddle point problem (1) of four different sizes are considered. In , it is stated that the iteration steps of both the SS and RSS preconditioned iteration methods remain constant
Figure 3. Eigenvalue distributions of the unpreconditioned saddle point matrix and of three SS-type preconditioned matrices ( ).
Table 3. Numerical results for three preconditioned GMRES methods with different grid sizes.
even if the problem size soars. Fortunately, this desirable property has also been carried over to GSS. Furthermore, GSS outperforms SS and RSS in that it requires the shortest CPU time and the least number of iteration steps in all test problems, as observed in Table 3. In light of this, we conclude that GSS is a competitive shift-splitting variant for solving the structured saddle point problem (1).
The shift-splitting iteration method/preconditioner for solving the block three- by-three saddle point problem (1) uses a single parameter to control the convergence. Being aware of the different roles of the block diagonal matrices in the saddle point matrix, we introduce an additional parameter to monitor the - block, which yields a generalized shift-splitting iteration method and a corresponding auxiliary preconditioner. By looking into the spectrum of the iteration matrix, we prove the unconditional convergence of the proposed method. We also apply the induced generalized shift-splitting preconditioner to speed up the convergence of GMRES. It is proved that the preconditioned matrix using the new preconditioner has a well-clustered spectrum which is verified by the numerical examples. Future works include finding the theoretically optimal choice of the parameters in the new scheme and develop other efficient shift-splitting preconditioners.
This work is supported by the National Natural Science Foundation (grant number: 11601323) and the Key Discipline Fund (2018) of College of Arts and Sciences in Shanghai Maritime University. The authors would like to thank the referees for their constructive suggestions.
1By “total iteration steps”, we mean the total number of matrix-vector products associated with ; for instance, if “iter”, as in Table 3, is denoted by 10(3), then the number of total iteration steps will be . In general, if “iter” is denoted by , then the number of total iteration steps will be .
 Bai, Z.-Z., Golub, G.H. and Ng, M.K. (2003) Hermitian and Skew-Hermitian Splitting Methods for Non-Hermitian Positive Definite Linear Sys-tems. SIAM Journal on Matrix Analysis and Applications, 24, 603-626.
 Bai, Z.-Z. (2015) Motiva-tions and Realizations of Krylov Subspace Methods for Large Sparse Linear Systems. Journal of Computational and Applied Mathematics, 283, 71-78. https://doi.org/10.1016/j.cam.2015.01.025
 Bramble, J.H., Pasciak, J.E. and Vassilev, A.T. (1997) Analysis of the Inexact Uzawa Algorithm for Saddle Point Prob-lems. SIAM Journal on Numerical Analysis, 34, 1072-1092. https://doi.org/10.1137/S0036142994273343
 Bai, Z.-Z. and Wang, Z.-Q. (2008) On Parameterized Inexact Uzawa Methods for Generalized Saddle Point Problems. Linear Algebra and Its Applications, 428, 2900- 2932. https://doi.org/10.1016/j.laa.2008.01.018
 Cao, Y. and Yi, S.-C. (2016) A Class of Uzawa-PSS Iteration Methods for Nonsingular and Singular Non-Hermitian Saddle Point Problems. Applied Mathematics and Computation, 275, 41-49. https://doi.org/10.1016/j.amc.2015.11.049
 Bai, Z.-Z. and Golub, G.H. (2007) Accelerated Hermitian and Skew-Hermitian Splitting Iteration Methods for Saddle-Point Problems. IMA Journal of Numerical Analysis, 27, 1-23. https://doi.org/10.1093/imanum/drl017
 Cao, Y., Yao, L., Jiang, M. and Niu, Q. (2013) A Relaxed HSS Preconditioner for Saddle Point Problems from Meshfree Discretization. Journal of Computational and Applied Mathematics, 31, 398-421. https://doi.org/10.4208/jcm.1304-m4209
 Benzi, M. and Guo, X.-P. (2011) A Dimensional Split Preconditioner for Stokes and Linearized Na-vier?-Stokes Equations. Applied Numerical Mathematics, 61, 66-76.
 Zhang, G.-F., Ren, Z.-R. and Zhou, Y.-Y. (2011) On HSS-Based Constraint Preconditioners for Generalized Saddle-Point Problems. Numerical Algorithms, 57, 273-287.
 Wang, N.-N. and Li, J.-C. (2019) A Class of New Extended Shift-Splitting Preconditioners for Saddle Point Problems. Journal of Computational and Applied Mathematics, 357, 123-145. https://doi.org/10.1016/j.cam.2019.02.015
 Cao, Y., Du, J. and Niu, Q. (2014) Shift-Splitting Preconditioners for Saddle Point Problems. Journal of Computational and Applied Mathematics, 272, 239-250.
 Salkuyeh, D.K., Masoudi, M. and Hezari, D. (2015) On the Generalized Shift-Splitting Preconditioner for Saddle Point Problems. Applied Mathematics Letters, 48, 55-61.
 Bai, Z.-Z. and Hadjidimos, A. (2014) Optimization of Extrapolated Cayley Transform with Non-Hermitian Positive Definite Matrix. Linear Algebra and Its Applications, 463, 322-339. https://doi.org/10.1016/j.laa.2014.08.021
 Zhang, K., Zhang, J.-L. and Gu, C.-Q. (2017) A New Relaxed PSS Preconditioner for Nonsymmetric Saddle Point Problems. Applied Mathematics and Computation, 308, 115-129. https://doi.org/10.1016/j.amc.2017.03.022
 Huang, N. and Ma, C.-F. (2019) Spectral Analysis of the Preconditioned System for the 3 × 3 Block Saddle Point Prob-lem. Numerical Algorithms, 81, 421-444.