Shift-splitting iteration methods for a class of large sparse linear matrix equations

1.   Introduction

Consider large sparse linear matrix equation

where $ A\in \mathbb{C}^{m \times m} $ and $ B\in \mathbb{C}^{n \times n} $ are non-Hermitian positive definite matrices, $ C\in \mathbb{C}^{m \times n} $ is a given complex matrix. In many areas of scientific computation and engineering applications, such as signal and image processing ^[9,20], control theory ^[11], photogrammetry ^[19], we need to solve such matrix equations. Therefore, solving such matrix equations by efficient methods is a very important topic.

We often rewrite the above matrix Eq (1.1) as the following linear system

where the vectors $ \mathbf{x} $ and $ \mathbf{c} $ contain the concatenated columns of the matrices $ X $ and $ C $, respectively, $ \otimes $ being the Kronecker product symbol and $ B^T $ representing the transpose of the matrix $ B $. Although this equivalent linear system can be applied in theoretical analysis, in fact, solving (1.2) is always costly and ill-conditioned.

So far there are many numerical methods to solve the matrix Eq (1.1). When the coefficient matrices are not large, we can use some direct algorithms, such as the QR-factorization-based algorithms ^[13,28]. Iterative methods are usually employed for large sparse matrix Eq (1.1), for instance, least-squares-based iteration methods ^[26] and gradient-based iteration methods ^[25]. Moreover, the nested splitting conjugate gradient (NSCG) iterative method, which was first proposed by Axelsson, Bai and Qiu in ^[1] to solve linear systems, was considered for the matrix Eq (1.1) in ^[14].

Bai, Golub and Ng originally established the efficient Hermitian and skew-Hermitian splitting (HSS) iterative method ^[5] for linear systems with non-Hermitian positive definite coefficient matrices. Subsequently, some HSS-based methods were further considered to improve its robustness for linear systems; see ^{[3,4,8,17,24,27]} and other literature. For solving the continuous Sylvester equation, Bai recently established the HSS iteration method ^[2]. Hereafter, some HSS-based methods were discussed for solving this Sylvester equation ^{[12,15,16,18,22,30,31,32,33]}. For the matrix Eq (1.1), Wang, Li and Dai recently use an inner-outer iteration strategy and then proposed an HSS iteration method ^[23]. According to the discussion in ^[23], if the quasi-optimal parameter is employed, the upper bound of the convergence rate is equal to that of the CG method. After that, Zhang, Yang and Wu considered a more efficient parameterized preconditioned HSS (PPHSS) iteration method ^[29] to further improve the efficiency for solving the matrix Eq (1.1), and Zhou, Wang and Zhou presented a modified HSS (MHSS) iteration method ^[34] for solving a class of complex matrix Eq (1.1).

Moreover, the shift-splitting (SS) iteration method ^[7] was first presented by Bai, Yin and Su to solve the ill-conditioned linear systems. Then this splitting method was subsequently considered for solving saddle point problems due to its promising performance; see ^[10,21] and other literature. In this paper, the SS technique is implemented to solve the matrix Eq (1.1). Some related convergence theorems of the SS method are discussed in detail. Numerical examples demonstrate that the SS is superior to the HSS and NSCG methods, especially when the coefficient matrices are ill-conditioned.

The content of this paper is arranged as follows. In Section 2 we establish the SS method for solving the matrix Eq (1.1), and then some related convergence properties are studied in Section 3. In Section 4, the effectiveness of our method is illustrated by two numerical examples. Finally, our brief conclusions are given in Section 5.

2.   The shift-splitting (SS) iteration method

Based on the shift-splitting proposed by Bai, Yin and Su in ^[7], we have the shift-splitting of $ A $ and $ B $ as follows:

and

where $ \alpha $ and $ \beta $ are given positive constants.

Therefore, using the splitting of the matrix $ A $ in (2.1), the following splitting iteration method to solve (1.1) can be defined:

Then, from the splitting of the matrix $ B $ in (2.2), we can solve each step of (2.3) iteratively by

Therefore, we can establish the following shift-splitting (SS) iteration method to solve (1.1).

Algorithm 1 (The SS iteration method). Given an initial guess $ X^{(0)}\in \mathbb{C}^{m \times n} $, for $ k = 0, 1, 2, \ldots $, until $ X^{(k)} $ converges.

Approximate the solution of

with $ R^{(k)} = C-AX^{(k)}B $, i.e., let $ Z^{(k)}: = Z^{(k, j+1)} $ and compute $ Z^{(k, j+1)} $ iteratively by

once the residual $ P^{(k)} = 2R^{(k)}-(\alpha I_m+A)Z^{(k, j+1)}B $ of the outer iteration (2.5) satisfies

where $ \|\cdot\|_F $ denotes the Frobenius norm of a matrix. Then compute

Here, $ \{\varepsilon_k\} $ is a given tolerance. In addition, we can choose efficient methods in the process of computing $ Z^{(k, j+1)} $ in (2.6).

The pseudo-code of this algorithm is shown as following:

Remark 1. Because the SS iteration scheme is only a single-step method, a considerable advantage is that it costs less computing workloads than the two-step iteration methods such as the HSS iteration ^[23] and the modified HSS (MHSS) iteration ^[34].

3.   Convergence analysis of the SS method

In this section, we denote by

the Hermitian and skew-Hermitian parts of the matrix $ A $, respectively. Moreover, $ \lambda_{\min} $ and $ \lambda_{\max} $ represent the smallest and the largest eigenvalues of $ H $, respectively, and $ \kappa = \lambda_{\max}/\lambda_{\min} $.

Firstly, the unconditional convergence property of the SS iteration (2.3) is given as follows.

Theorem 1. Let $ A\in \mathbb{C}^{m \times m} $ be positive definite, and $ \alpha $ be a positive constant. Denote by

Then the convergence factor of the SS iteration method (2.3) is given by the spectral radius $ \rho(M(\alpha)) $ of the matrix $ M(\alpha) $, which is bounded by

Consequently, we have

i.e., the SS iteration (2.3) is unconditionally convergent to the exact solution $ X^{\star}\in \mathbb{C}^{m \times n} $ of the matrix Eq (1.1).

Proof. The SS iteration (2.3) can be reformulated as

Using the Kronecker product, we can rewrite (3.3) as follows:

where $ M(\alpha) $ is the iteration matrix defined in (3.1), and $ N(\alpha) = 2B^{-T}\otimes(\alpha I_m+A)^{-1} $.

We can easily see that $ \rho(M(\alpha))\leq\varphi(\alpha) $ holds for all $ \alpha > 0 $. From Lemma 2.1 in ^[8], we can obtain that $ \varphi(\alpha) < 1 $, $ \forall \alpha > 0 $. This completes the proof.

Noting that the matrix $ (\alpha I_m+A)^{-1}(\alpha I_m-A) $ is an extrapolated Cayley transform of $ A $, from ^[6], we can obtain another upper bound for the convergence factor of $ \rho(M(\alpha)) $, as well as the minimum point and the corresponding minimal value of this upper bound.

Theorem 2. Let the conditions of Theorem 1 be satisfied. Denote by

Then for the convergence factor of $ \rho(M(\alpha)) $ it holds that

Moreover, at

the function $ \phi(\alpha) $ attains its minimum

where

Proof. From Theorem 3.1 in ^[6], we can directly obtain (3.2)–(3.7).

Remark 2. $ \alpha_\star $ is called the theoretical quasi-optimal parameter of the SS iteration method. Similarly, the theoretical quasi-optimal parameter $ \beta_\star $ of the inner iterations (2.6) can be also obtained, which has the same form as $ \alpha_\star $.

In the following, we present another convergence theorem for a new form.

Theorem 3. Let the conditions of Theorem 1 be satisfied. If $ \{X^{(k)}\}_{k = 0}^{\infty}\subseteq \mathbb{C}^{m \times n} $ is an iteration sequence generated by Algorithm 1 and if $ X^{\star}\in \mathbb{C}^{m \times n} $ is the exact solution of the matrix Eq (1.1), then it holds that

where the constants $ \mu $ and $ \theta $ are given by

In particular, when

the iteration sequence $ \{X^{(k)}\}_{k = 0}^{\infty} $ converges to $ X^{\star} $, where $ \varepsilon_{\max} = \max_{k}\{\varepsilon_k\} $.

Proof. We can rewrite the SS iteration in Algorithm 1 as the following form:

with $ \mathbf{r}^{(k)} = \mathbf{c}-(B^{T}\otimes A)\mathbf{x}^{(k)} $, where $ \mathbf{z}^{(k)} $ is such that the residual

satisfies $ \|\mathbf{p}^{(k)}\|_2\leq\varepsilon_k\|\mathbf{r}^{(k)}\|_2 $.

In fact, the inexact variant of the SS iteration method for solving the linear system (1.2) is just the above iteration scheme (3.9). From (3.9), we obtain

Because $ \mathbf{x}^{\star}\in \mathbb{C}^{n} $ is the exact solution of the linear system (1.2), it must satisfy

By subtracting (3.11) from (3.10), we have

Taking norms on both sides from (3.12), then

Noticing that

by (3.13) the estimate

can be obtained. Note that for a matrix $ Y\in \mathbb{C}^{m \times n} $, $ \|Y\|_F = ||\mathbf{y}||_2 $, where the vector $ \mathbf{y} $ contains the concatenated columns of the matrix $ Y $. Then the estimate (3.14) can be equivalently rewritten as

So we can easily get the above conclusion.

Remark 3. From Theorem 3 we know that, in order to guarantee the convergence of the SS iteration, it is not necessary for the condition $ \varepsilon_k\rightarrow 0 $. All we need is that the condition (3.8) is satisfied.

4.   Numerical results

In this section, two different matrix equations are solved by the HSS, SS and NSCG iteration methods. The efficiencies of the above iteration methods are examined by comparing the number of outer iteration steps (denoted by IT-out), the average number of inner iteration steps (denoted by IT-in-1 and IT-in-2 for the HSS, IT-in for the SS), and the elapsed CPU times (denoted by CPU). The notation "–" shows that no solution has been obtained after 1000 outer iteration steps.

The initial guess is the zero matrix. All iterations are terminated once $ X^{(k)} $ satisfies

We set $ \varepsilon_k = 0.01 $, $ k = 0, 1, 2, \ldots $ to be the tolerances for all the inner iteration schemes.

Moreover, in practical computation, we choose direct algorithms to solve all sub-equations involved in each step. We use Cholesky and LU factorization for the Hermitian and non-Hermitian coefficient matrices, respectively.

Example 1 (^[2]) We consider the matrix Eq (1.1) with $ m = n $ and

where $ M, N\in \mathbb{R}^{n \times n} $ are two tridiagonal matrices as follows:

In Tables 1 and 2, the theoretical quasi-optimal parameters and experimental optimal parameters of HSS and SS are listed, respectively. In Tables 3 and 4, the numerical results of HSS and SS are listed.

From Tables 3 and 4 it can be observed that, the SS outperforms the HSS for various $ n $ and $ q $, especially when $ q $ is small (the coefficient matrices are ill-conditioned).

Moreover, as two single-step methods, the numerical results of NSCG and SS are compared in Table 5. From Table 5 we see that the SS method has better computing efficiency than the NSCG method.

Example 2 (^[2]) We consider the matrix Eq (1.1) with $ m = n $ and

where $ L $ is a strictly lower triangular matrix and all the elements in the lower triangle part are ones, and $ t $ is a specified problem parameter. In our tests, we take $ t = 1 $.

In Tables 6 and 7, for various $ n $ and $ r $, we list the theoretical quasi-optimal parameters and experimental optimal parameters of HSS and SS, respectively. In Tables 8 and 9, the numerical results of HSS and SS are listed. Moreover, the numerical results of NSCG and SS are compared in Table 10.

From Tables 8–10 we get the same conclusion as example 1.

Therefore, for large sparse matrix equation $ AXB = C $, the SS method is an effective iterative approach.

5.   Conclusions

By utilizing an inner-outer iteration strategy, we established a shift-splitting (SS) iteration method for large sparse linear matrix equations $ AXB = C $. Two different convergence theories were analysed in depth. Furthermore, the quasi-optimal parameters of SS iteration matrix are given. Numerical experiments illustrated that, the SS method can always outperform the HSS and NSCG methods both in outer and inner iteration numbers and computing time, especially for the ill-conditioned coefficient matrices.

Acknowledgments

The authors are very grateful to the anonymous referees for their helpful comments and suggestions on the manuscript. This research is supported by the Natural Science Foundation of Gansu Province (No. 20JR5RA464), the National Natural Science Foundation of China (No. 11501272), and the China Postdoctoral Science Foundation funded project (No. 2016M592858).

Conflict of interest

The authors declare there is no conflict of interest.