An efficient relaxed shift-splitting preconditioner for a class of complex symmetric indefinite linear systems

1.   Introduction

Consider the iterative solution of the following complex linear equations of the form

where $\widetilde{\mathscr{A}} = W+iT\in \mathbb{C}^{n\times n}$ is a complex matrix, with $W, \; T \in \mathbb{R}^{n\times n}$ being symmetric matrices, $\tilde{x} = y+iz, \; \tilde{b} = f+ig$, and $i = \sqrt{-1}$ denotes the imaginary unit.

Complex systems such as (1.1) are important and arise in various scientific computing and engineering applications, such as diffuse optical tomography, structural dynamics ^[7], optimal control problems for PDEs with various kinds of state stationary or time dependent equations, e.g., Poisson, convection diffusion, Stokes ^[2], wave propagation and so on. More details on this class of questions are given in references ^{[1,3,5,14,15]}.

When the matrices $W, \; T$ are symmetric positive semi-definite with at least one of them being positive definite, Bai et al.^[7,8] proposed the modified Hermitian and skew-Hermitian splitting (MHSS) iteration method and the preconditioned MHSS (PMHSS) iteration methods to compute an approximate solution for the complex linear systems (1.1), see also ^[6]. It is proved in ^[7,8] that the MHSS and PMHSS methods converge to the unique solution of (1.1) unconditionally. Moreover, Bai et al. pointed out the $h$-independent behavior of the corresponding preconditioner of the PMHSS iteration method. To solve the systems (1.1) further and more efficiently, Zheng et al. ^[24] designed a double scale splitting (DSS) iteration method and also analyzed the unconditional convergence property. Furthermore, two reciprocal optimal iteration parameters and corresponding optimal convergence factor are determined simultaneously. There are some other effective iteration methods at the same time, such as Euler preconditioned SHSS iteration method ^[17], Double parameter splitting (DPS) iteration method ^[18], etc.

However, the matrix $T$ of complex symmetric system of linear equations arises in direct frequency domain analysis ^[10] and time integration of parabolic partial differential equations ^[4] is usually symmetric indefinite, the MHSS, PMHSS and DSS methods may be applicative or not, due to the fact that the coefficient matrices $\alpha I+T$, $\alpha V+T$ and $\alpha W+T$ are indefinite or singular. For such a problem, multiply the complex linear systems on the left by $T$, Wu ^[19] developed the simplified Hermitian normal splitting (SHNS) iteration method. In order to accelerate the convergence of the SHNS method, Zhang et al. ^[21] established a preconditioned SHNS (PSHNS) iteration method and constructed a corresponding preconditioner. Although these two iteration methods are unconditionally convergent, they still involve the complex arithmetics in each inner iteration, which can result in expensive computational costs. More importantly, computation of the optimal values with any of the aforementioned two methods is a time-consuming because it first needs to compute the maximum and the minimum eigenvalues of some dense matrices.

In this paper, we will focus on the case that $W$ is symmetric positive definite and $T$ is symmetric indefinite. In order to avoid the complex arithmetic, the complex linear systems (1.1) are often transformed into the real block two-by-two systems as follows ^[3,25]:

this real form can be regarded as a special class of generalized saddle point problems ^[11]. Based on the relaxed preconditioning technique ^[12] for generalized saddle point problems, Zhang et al. ^[22] proposed a block splitting (BS) preconditioner. In order to overcome the nonzero off-diagonal block becoming unbounded as the relaxed parameter approaches 0, Zhang et al. ^[23] proposed an improved block (IB) splitting preconditioner. All of these preconditioners are highly close to the coefficient matrix of the real linear systems (1.2), when accelerating the Krylov subspace methods, a linear sub-system with coefficient matrix $\alpha W+T^2$ must be solved in each inner iteration. However, $\alpha W+T^2$ is a dense symmetric positive and definite matrix. Unlike sparse matrices, the computation of a dense matrix is more difficult, for some high dimensional problems, it may be hard to solve.

Fortunately, constructing these preconditioners give us some inspiration, while the linear systems (1.2) need to make some changes first. We introduce a scalar matrix $\alpha I$ to construct a new complex block two-by-two linear systems, then based on the shift-splitting preconditioner ^[13] for saddle point problems, we propose an efficient relaxed shift-splitting (ERSS) preconditioner. This preconditioner not only avoid the complex arithmetics but also maintain the sparse properties of the matrices $W$ and $T$. More importantly, the ERSS preconditioner is highly close to the original coefficient matrix of the new complex block two-by-two linear systems and the relax parameter is easily to implemented.

The remainder of this work is organized as follows. In Section 2.1, we propose an efficient relaxed shift-splitting (ERSS) preconditioner and the eigenvalue properties of preconditioned matrix are discussed. In Section 2.2, using the scaled norm minimization (SNM) method ^[20], we derive a practical formula for computing the parameter value $\alpha$. In Section 3, numerical experiments are presented to show the effectiveness of the ERSS preconditioner. Finally, we end this paper with some conclusions in Section 4.

2.   The ERSS preconditioner

In this section, we first build the ERSS preconditioner and derive some spectral properties of the corresponding preconditioned system. Then by using the scaled norm minimization (SNM) method, a practical estimation formula is given to compute the relaxed parameter.

2.1.   The ERSS preconditioner

Firstly, we reconstruct the complex linear systems (1.1) into the following structure by introducing a scalar matrix $\alpha I$:

where $\alpha$ is a positive constant. We regard this system (2.1) as a "saddle point system". Based on the shift-splitting preconditioner ^[13] for saddle point problems, we propose the following relax shift-splitting preconditioner:

The difference between $\mathscr{P}_{\rm{ERSS}}$ and $\mathscr{A}$ is

Only the (1, 2) block being nonzero in $\mathscr{R}_{\rm{ERSS}}$ shows that the preconditioner $\mathscr{P}_{\rm{ERSS}}$ is a good approximation to the coefficient matrix $\mathscr{A}$ and it may be easier to analyze the eigenvalue distributions of the preconditioned matrix $\mathscr{P}_{\rm{ERSS}}^{-1}\mathscr{A}$.

In actual implementations, the actions of the preconditioned Krylov subspace methods with the preconditioner $\mathscr{P}_{\rm{ERSS}}$, are often realized through solving a sequence of generalized residual equations of the form $\mathscr{P}_{\rm{ERSS}}z = r$, where $r = (r_1^*, r_2^*)^*\in\mathbb{C}^{2n}$, with $r_1, \; r_2\in\mathbb{C}^n$ represent the current residual vector, $z = (z_1^*, z_2^*)^*\in\mathbb{C}^{2n}$, with $z_1, \; z_2\in\mathbb{C}^n$ represent the generalized residual vector, i.e.,

By using the matrix factorization of $\mathscr{P}_{\rm{ERSS}}^{-1}$, we obtain the following procedure for the residual vector $z = (z_1^*, \; z_2^*)^*$:

Algorithm 1:

(1) solve $(\alpha I+\frac{1}{\alpha}W)u_1 = r_2-\frac{1}{\alpha}W r_1;$

(2) $z_1 = \frac{1}{\alpha}(r_1+u_1);$

(3) solve $Tu_2 = u_1;$

(4) $z_2 = -i\alpha u_2.$

From Algorithm 1, we see that two linear subsystems with sparse real coefficient matrices $\alpha I+\frac{1}{\alpha}W$ and $T$ need to be solved at steps (1) and (3). Since the matrix $\alpha I+\frac{1}{\alpha}W$ is symmetric positive and definite and $T$ is symmetric, both of them are sparse matrices, then the above linear subsystems can be solved effectively by sparse Cholesky factorization and LU method, respectively.

The spectral distributions of the preconditioned matrix relate closely to the convergence rate of Krylov subspace methods ^[9]. The following result shows the eigenvalue distributions of the preconditioned matrix $\mathscr{P}_{\rm{ERSS}}^{-1}\mathscr{A}$.

Theorem 2.1. Let $W\in\mathbb{R}^{n\times n}$ is symmetric and positive definite and $T\in\mathbb{R}^{n\times n}$ is symmetric indefinite, $\alpha$ is a positive constant. Then the preconditioned matrix $\mathscr{P}_{\mathit{{ERSS}}}^{-1}\mathscr{A}$ has eigenvalues at 1 with multiplicity $n$, and the remaining $n$ eigenvalues are of the form $\frac{\alpha^2 (\omega -i\tau)}{1+\alpha^2 \omega}$, where

Proof. The preconditioned matrix can be rewritten as

From (2.2), we can get

where $\Theta = i T^{-1}(\alpha I+\frac{1}{\alpha}W)^{-1} W (\alpha I-\frac{i}{\alpha}T)$. Define $\widetilde{\Theta} = (\frac{1}{\alpha}I+\alpha W^{-1})^{-1}(\frac{1}{\alpha}I+i\alpha T^{-1})$, then $\Theta$ is similar to $\widetilde{\Theta}$. Assume that $(\tilde{\lambda}, \; \mu)$ is an eigenpair of $\widetilde{\Theta}$, i.e., $\widetilde{\Theta}\mu = \tilde{\lambda} \mu$,

Multiplying $\frac{\mu^*}{\mu^*\mu}$ by the right and left side of Eq (2.3), by simple calculations, we have

Hence the eigenvalues of the preconditioned matrix $\mathscr{P}_{\text{ERSS}}^{-1}\mathscr{A}$ are at 1 with multiplicity $n$, and the remaining $n$ eigenvalues are of the form $\frac{\alpha^2 (\omega -i\tau)}{1+\alpha^2 \omega}$.

Remark 2.1. Let $W\in\mathbb{R}^{n\times n}$ is symmetric and positive definite and $T\in\mathbb{R}^{n\times n}$ is symmetric indefinite, $\alpha$ is a positive constant. Then the non-unit eigenvalues of the preconditioned matrix $\mathscr{P}_{\text{ERSS}}^{-1}\mathscr{A}$ are clustered at $0_+$ if $\alpha$ is close to 0, while the real parts of the eigenvalues are around at $1$ if $\alpha$ approaches to $\infty$.

Remark 2.2. If $|\tau|_{\max}\leq \omega_{\min}$, then for $\forall \alpha > 0$, all eigenvalues of $\mathscr{P}_{\text{ERSS}}^{-1}\mathscr{A}$ satisfies $|\lambda-1| < 1$, where $|\tau|_{\max}$ is the maximum value of the absolute value of eigenvalues of the matrix $T^{-1}$ and $\omega_{\min}$ is the smallest eigenvalue of the matrix $W^{-1}$.

2.2.   Practical estimation for the parameter $\alpha$

When $\mathscr{P}_{\text{ERSS}}$ is used as a preconditioner, we expect that $\mathscr{P}_{\text{ERSS}}$ is as close as possible to the coefficient matrix $\mathscr{A}$ of the complex linear systems (2.1). So we try to derive a practical formula for computing the optimal parameter $\alpha$ such that $\mathscr{R}_{\text{ERSS}}\approx 0$. Recently, Yang ^[20] proposed an easily implemented scaled norm minimization (SNM) method to compute the parameter values including several traces of some matrices for the Hermitian and skew-Hermitian splitting (HSS) method ^[9,14,16]. Here, we define $tr(\cdot)$ as a matrix's trace. Owing to $\|\mathscr{R}_{\text{ERSS}} \|_{F}^{2} = tr(\mathscr{R}_{\text{ERSS}}^*\mathscr{R}_{\text{ERSS}})$, we first give

Because $tr(A+B) = tr(A)+tr(B)$ and $tr(k\cdot A) = k\cdot tr(A)$ for any $A\in\mathbb{R}^{n\times n}$ and $k\in\mathbb{R}$. It follows that

It's clear to know when $\alpha_\star = \frac{\sqrt{\|T\|_F}}{\sqrt[4]{n}}$ that minimizes $\|\mathscr{R}_{\text{ERSS}}\|_{F}^{2}$. Obviously, the calculation of relax parameter $\alpha_\star$ is easy to realise.

3.   Numerical examples

In this section, we employ two examples to test the performances of the ERSS preconditioner in terms of both iteration count (denoted as IT) and computing time (in seconds, denoted as CPU). To show the effectiveness of the ERSS preconditioner (2.2), we also test the other two preconditioners: PSHNS preconditioner $\mathscr{P}_{\text{PSHNS}}$ ^[21] and IB preconditioner $\mathscr{P}_{\text{IB}}$ ^[23] as follows:

$\mathscr{P}_{\text{PSHNS}}$ is used to precondition the complex linear systems (1.1), and $\mathscr{P}_{\text{IB}}$ preconditions the real linear systems (1.2).

In implementations, we use those preconditioners to accelerate the convergence of the generalized minimum residual (GMRES) method. The initial guess $x^{(0)}$ for the preconditioned GMRES method is chosen to zero vector, and the iterations are terminated once the current iterate $x^{(k)}$ satisfies

The relaxed parameters used in both PSHNS and IB preconditioners are the experimentally found ones, which minimize the number of iteration steps, while the relaxed parameter of ERSS preconditioner is $\alpha_\star = \sqrt{\|T\|_F}/\sqrt[4]{n}$. In addition, the systems of linear equations involved in the preconditioned GMRES method are solved by direct methods, that is, the Cholesky factorization in combination with the symmetric approximate minimum degree reordering and the LU factorization in combination with the column approximate minimum degree reordering, respectively.

All experiments are performed by using MATLAB (version R2009b) in double precision on a personal computer with 3.60GHz central processing unit (Intel(R) Core(TM) i7-4790 CPU), 8.00G memory and Windows 7 operating system.

Example 3.1 (See ^[8,23,24,25]) The following complex symmetric linear system is considered

where $M$ and $K$ are the inertia and stiffness matrices, respectively; $C_V$ and $C_H$ are the viscous and hysteretic damping matrices, respectively; and $\omega$ is the driving circular frequency.

In our numerical computations, we take $C_H = 0.02 K$, $\omega = 2\pi$, $C_V = \frac{1}{2}M$, $M = kI$ and $K$ is the five-point centered difference matrix approximating the negative Laplacian operator with homogeneous Dirichlet boundary conditions, on a uniform mesh in the unit square $[0, \; 1]\times [0, \; 1]$ with the mesh size $h = \frac{1}{m+1}$. In this case, the matrix $K\in\mathbb{R}^{n\times n}$ possesses the tensor-product form $K = I\otimes V_m+ V_m\otimes I$ with $V_m = h^{-2}\text{tridiag}(-1, 2, -1)\in\mathbb{R}^{m\times m}$. Hence, the total number of variables is $n = m^2$. In addition, the right-hand side vector $\tilde{b} = (1+i)\widetilde{\mathscr{A}}*\text{ones}$(n, 1). Furthermore, we normalize the coefficient matrix and right-hand side by multiplying both by $h^2$.

In Table 1, we report results for GMRES preconditioned with ERSS, PSHNS and IB preconditioners for different mesh-size $h$ and symmetric positive and definite matrices $M$. From these results we observe that when used as a preconditioner, by choosing the theoretical optimal parameter $\alpha_\star$, ERSS performs much better than PSHNS and IB in both iteration steps and CPU times, especially when the mesh-size $h$ becomes small. While the number of iterations with the PSHNS and IB preconditioners increase with problem size, those for the ERSS preconditioner are almost constant. In addition, searching for optimal parameters of the PSHNS and IB preconditioners is quite time-consuming, especially for the latter, while the calculation of the parameter of the ERSS preconditioner is effortless.

Example 3.2. (See ^[23]) Consider the system of linear Eq (1.1) as following:

where $K = I_m\otimes V_m+V_m\otimes I_m$, $\tau = 2\pi^2$, $h = \frac{1}{m+1}$, $n = m^2$ and $V_m = h^{-2}tridiag(-1, 2, -1)\in\mathbb{R}^{m\times m}$ is a tridiagonal matrix. $\omega = \sqrt{k}\pi^2$ is a variable.

We choose the symmetric positive and definite matrix $W = K+(3+\sqrt{3})\tau I_{m^2}$ and the symmetric indefinite matrix $T = K-(3-\sqrt{3})\omega I_{m^2}$, the right-hand side vector $\tilde{b} = (1+i)\widetilde{\mathscr{A}}*\text{ones}(m^2, 1)$. Furthermore, we normalize the coefficient matrix and right-hand side by multiplying both by $h^2$.

In Table 2, we list results for GMRES preconditioned with ERSS, PSHNS and IB for different mesh-size $h$ and variable $k$. Note that the parameter values of the PSHNS and IB preconditioners are the experimentally found optimal ones, which is time-consuming. Despite the iteration steps and CPU times of ERSS preconditioner exceed that of IB as the mesh-size $h$ decreases, the difference is acceptable. As a consequence, although the number of iteration steps of the ERSS method increase slightly with the mesh refinement, it still has strong competitiveness due to the fast parameter calculation method.

Finally, we present the experimental optimal results for the ERSS-preconditioned GMRES method by minimizing the numbers of iterations with respect to different test examples and variables in Table 3. We can see that the experimental optimal parameters in Table 3 are consistent with theoretical optimal parameters $\alpha_\star = \sqrt{\|T\|_F}/\sqrt[4]{n}$ in Tables 1 and 2. While, the experimental optimal results are slightly larger than that of theoretical optimal results in Table 2, and this insignificant difference in iteration steps are acceptable. Table 3 demonstrates that the ERSS-preconditioned GMRES method is efficient and stable when the relaxation parameter selected theoretically optimal $\alpha_\star$.

Eigenvalue distributions ($48\times 48$ grids) of the three preconditioned matrices are plotted in Figures 1 and 2 for different variables $k$. It is evident that the ERSS preconditioned matrix is of a well-clustered spectrum around $1$ away from zero, especially in Example 3.1.

4.   Conclusions

To solve a class of complex linear systems (1.1), an efficient relaxed shift-splitting preconditioner is proposed in this paper by introducing a scalar matrix $\alpha I$. The new preconditioner not only remains easy computational but also is closer to the original two-by-two complex coefficient matrix (2.1). Theoretical analysis proves that the preconditioned matrix has a well-clustered eigenvalue distribution with a reasonable choice of the relaxation parameters. More importantly, an efficient and practical formula for computing the relax parameter value $\alpha$ is derived by computing dimension and Frobenius norm of the matrix $T$. Numerical experiments are presented to illustrate that the presented preconditioner is feasible and effective compared with other existing block preconditioners.

Acknowledgments

This work is supported by the Teacher education curriculum reform research project in Henan Province (NO.2021JSJYYB124), Xinyang City Federation of Social Science planning project (NO.2021JY063), Xinyang University scientific research project (NO.2022-XJLYB-002).

Conflict of interest

The authors declare no conflict of interest.