Two accelerated gradient-based iteration methods for solving the Sylvester matrix equation AX + XB = C

1.   Introduction

In this paper, we consider the iterative solution of the following Sylvester matrix equation:

where $A\in \mathbb{R}^{m\times m}$, $B\in \mathbb{R}^{n\times n}$, $C\in \mathbb{R}^{m\times n}$ are constant matrices and $X\in \mathbb{R}^{m\times n}$ is the unknown matrix to be obtained.

Due to the extensive applications of Eq (1.1) in control theory and stability analysis ^[1,10,14], it has garnered considerable attention, and many algorithms have been proposed over the past few decades. For example, the gradient-based iteration (GI) algorithm described in ^[6,8,9] has proven to be an effective method for solving Eq (1.1). By incorporating a tunable parameter into the GI algorithm, a relaxed gradient-based iteration (RGI) algorithm ^[18] was introduced, which demonstrates improved performance over the GI algorithm when the relaxed parameter is appropriately adopted. To enhance the RGI algorithm's efficiency, an accelerated gradient-based iteration (AGBI) algorithm was proposed by leveraging the latest information from the preceding half-step in ^[25]. In order to achieve a lower computational cost, a Jacobi gradient iteration (JGI) method was outlined in ^[13] based on the Jacobi splitting of $A$ and $B$. Drawing inspiration from the AGBI and JGI algorithms, Tian et al. ^[21] further developed an accelerated JGI (AJGI) algorithm. Additionally, various other iteration algorithms ^[7,20,22] have been devised for solving Eq (1.1) and other related matrix equations ^{[17,23,24,29]}, because of its wide applications.

Preconditioning techniques aim to alter the spectral characteristics of matrices through linear transformations, which are often integrated with other iteration methods and lead to various new algorithms such as the preconditioned HSS method ^[2,16], generalized preconditioned HSS methods ^[28], and preconditioned MHSS iteration methods ^[4], etc. The heavy-ball momentum method is widely applied to accelerate the convergence rate of the gradient method ^[5,19]. In this paper, inspired by the references ^[2,5,19], we combine the precondition technique and the momentum item with the gradient-based iteration algorithm, and the specific work can be summarized as follows:

(a) Novel Methodology. We have developed the preconditioned gradient-based iteration (PGI) and gradient-based momentum iteration (GMI) algorithms for solving Eq (1.1), which are more efficient than existing methods in terms of computational complexity and accuracy.

(b) Theoretical Insights. Our work provides new theoretical insights into gradient-based iteration algorithms. The convergences of PGI and GMI algorithms are rigorously proved.

(c) Adaptive Parameter Selection. We have developed a new parameter selection strategy that minimizes the current residual norm, leading to improved performance of the proposed algorithms. This strategy is practical and can be easily implemented in various numerical algorithms, enhancing their efficiency and accuracy.

(d) Empirical Results. Through extensive numerical experiments, we have shown that our methods outperform current state-of-the-art techniques in the solving of Eq (1.1).

The remainder of this paper is organized as follows: In Section 2, we first review the GI algorithm, and present the PGI and GMI algorithms, whose convergence properties are analyzed in detail. In Section 3, we construct the adaptive PGI and GMI algorithms in which the parameters are updated by utilizing the iterative information. In Section 4, several numerical examples are employed to show the robustness and efficiencies of the proposed algorithms. Finally, some conclusions are drawn in the last section.

2.   Two accelerated GI algorithms

In this section, we first review the GI algorithm. Subsequently, two accelerated GI algorithms are presented, and detailed analyses are conducted on their convergence properties. In the following, several lemmas are given, which will be used in the subsequent proofs.

Lemma 2.1. ^[12] Let $A\in \mathbb{R}^{m\times n}, B\in \mathbb{R}^{p\times q}$, and

For any matrix $M\in \mathcal{R}(A, B)$, it holds that

where $\sigma_{min}(A)$ and $\sigma_{min}(B)$ are the smallest singular values of the matrices $A$ and $B$, respectively.

Lemma 2.2. ^[11,15] Both roots of the real quadratic equation $x^2-bx+c = 0$ are less than one in modulus if and only if $|c| < 1$ and $|b| < 1+c$.

By utilizing the hierarchical identification principle, Eq (1.1) can be reformulated into two subsystems as follows:

The GI algorithm for solving (2.1) can be described as follows:

It is shown that the GI algorithm ^[8] converges when

where $\lambda_{max}(AA^T)$ and $\lambda_{max}(B^TB)$ are the largest eigenvalues of $AA^T$ and $B^TB$, respectively.

2.1.   The PGI algorithm

By introducing two preconditioners, $P$ and $Q$, in Algorithm 1, a preconditioned gradient-based iterative (i.e., PGI) algorithm is constructed and summarized as follows:

Remark 1. If $P = I_{m}$ and $Q = I_{n}$ are adopted, the PGI iteration method is reduced to the original GI algorithm in ^[8], where $I_s$ is an identity matrix with size $s$.

Remark 2. Two practical choices of the matrices $P$ and $Q$ are listed as follows:

1) $P = {\rm diag}(A), Q = {\rm diag}(B),$ where ${\rm diag}(A)$ and ${\rm diag}(B)$ are the diagonal matrices of $A$ and $B$, respectively.

2) $P = {\rm tridiag}(A^{T}A), Q = {\rm tridiag}(BB^{T}),$ where ${\rm tridiag}(A^{T}A)$ and ${\rm tridiag}(BB^{T})$ are the tridiagonal matrices of $A^{T}A$ and $BB^{T}$, respectively.

Theorem 2.1. Let $X^*$ be the solution of Eq (1.1). The iterative solution $X^{(k)}$ generated by Algorithm 2 converges to $X^*$ for any initial value if and only if the parameter $\mu$ satisfies the condition

where $\|\cdot\|_2$ is the 2-norm of the matrix.

Proof: For $k = 1, 2, \cdots$, define the $k$th error matrices $\widetilde{X}^{(k)}: = X^{(k)}-X^*$, which satisfy the following recurrence:

By using the Kronecker product ^[26], the above equation can be reformulated by

Taking the 2-norm of $vec(\widetilde{X}^{(k)})$, it follows that

where

Thus,

If $\mu$ satisfies (2.2), we know that $\widetilde{X}^{(k)}\rightarrow 0$ as $k\rightarrow \infty$. The proof is completed.

2.2.   The GMI algorithm

In order to make full use of the information from the previous iteration step, a momentum term will be added to the GI algorithm, and then the second accelerated GI (i.e., GMI) algorithm will be proposed and summarized as follows:

Remark 3. If $\beta$ is chosen to be 0, then the GMI algorithm is just the GI algorithm.

Theorem 2.2. Assume that the matrices $A$ and $B$ are non-singular. If Eq (1.1) has a unique solution $X^*$, then the iterative solution $X^{(k)}$ obtained from Algorithm 3 converges to $X^*$ for any initial values if and only if the parameters $\mu$ and $\beta$ satisfy the following conditions:

Case 1: When $\ 3q_2-q_1^2 > 0$,

or

where $q_1 = \sigma_{min}^2(A)+\sigma_{min}^2(B)-2\|A\|_{2}\|B\|_{2}$, $q_2 = (\|A\|_{2}^2+\|B\|_{2}^2)(\|A\|_{2}+\|B\|_{2})^2$, $a = \min\left\{\frac{1}{2}, \sqrt[]{\frac{q_1^2-q_2}{8q_2}}\right\}$, $b = \min\left\{\frac{1}{2}, \sqrt{\frac{q_1^2-2q_2}{4q_2}}\right\}$, $c = q_1^2-q_2(8\beta^2+1)$ and $d = q_1^2-q_2(4\beta^2+2)$.

Case 2: When $\ 3q_2-q_1^2\leq0$,

Proof: Define the error matrices

From Algorithm 3, it follows that

Taking the $F$-norm on (2.7), it yields that

and

By the triangle inequality and the property of the $F$-norm, we have

or

Squaring both sides, we obtain:

According to Lemma 2.1, we have

Combining (2.8)–(2.10), it yields that

By (2.11), we have

where

Let $\lambda$ be the eigenvalue of the matrix $H$. We know that

It then follows from Lemma 2.2 that $|\lambda| < 1$ if and only if

which implies that

In addition, $H\geq0$ ($H\geq0$, if $h_{ij}\geq0$ holds for all $1 < i < 2$, $1 < j < 2$.) if and only if

Together with (2.12) and (2.13), we have

In the following, we mainly solve the inequalities (2.14) to obtain the range of $\mu$ and $\beta$. The second inequality of (2.14) is equivalent to

Consider (2.15) as a system of quadratic inequalities in terms of $\mu$, and determine the range of $\mu$ to make the inequalities hold. Let's first solve the first inequality of (2.15). When $\Delta_1 = 4q_1^2-4q_2(8\beta^2+1) > 0$, i.e., $0 < \beta < \sqrt[]{\frac{q_1^2-q_2}{8q_2}},$ there are two solutions $\frac{q_1-\sqrt{q_1^2-q_2(8\beta^2+1)}}{q_2}$ and $\frac{q_1+\sqrt{q_1^2-q_2(8\beta^2+1)}}{q_2}$ for the quadratic equation $q_2\mu^2-2q_1\mu +8\beta^2+1 = 0$. So the solution of the first inequality in (2.15) is

When $\Delta_1\leq0$, $q_2\mu^2-2q_1\mu +8\beta^2+1$ is always greater than or equal to 0. So the inequality of $q_2\mu^2-2q_1\mu +8\beta^2+1 < 0$ has no solution.

Solving the the second inequality of (2.15) by the same method in the following. When $\Delta_2 = 4q_1^2-4q_2(4\beta^2+2) \leq 0$, i.e., $\beta\geq\sqrt{\frac{q_1^2-2q_2}{4q_2}},$ $q_2\mu^2-2q_1\mu +4\beta^2+2$ is always greater than or equal to 0. So the solution of the second inequality of (2.15) is $\mu\in R;$ when $\Delta_2 > 0$, i.e., $0 < \beta < \sqrt{\frac{q_1^2-2q_2}{4q_2}},$ the solution of the second inequality of (2.15) is

In order to find the solution for (2.15), we need to consider the following cases:

Case 1: If $3q_2-q_1^2 > 0$, then

or

Case 2: If $3q_2-q_1^2\leq0$, then

Together with the first inequality of (2.14) and (2.16)–(2.18), (2.4)–(2.6) are obtained. Thus, the proof is completed. □

Remark 4. When the error iteration matrix of the proposed algorithm is of size $2 \times 2$, the idea of using Lemma 2.2 to determine the range of the coefficients in quadratic equations, thereby ensuring the convergence of the algorithm, has been widely applied in many literatures. For example, the SOR-like methods for solving the absolute value equations ^[11,15].

3.   The APGI and AGMI algorithms

In this section, explicitly giving the varied parameters in the proposed algorithms by minimizing the residual in every iteration, the PGI and GMI algorithms with adaptive parameters are constructed.

3.1.   The adaptive PGI algorithm

We first present the calculation rule for the parameter used in Algorithm 2 by minimizing the current residual norm. The details are described below:

Suppose the parameter $\mu_k$ is taken in Algorithm 2 and for $k = 1, 2, 3, \cdots$, the previous $k-1$ residuals are defined by

then the PGI algorithm can be simply rewritten as

According to (3.1) and (3.2), the $k$th residual is given by

with $M^{(k-1)} = P^{-1}A^TR^{(k-1)}B+R^{(k-1)}B^TQ^{-1}B+AP^{-1}A^TR^{(k-1)}+AR^{(k-1)}B^TQ^{-1}.$

Taking the $F$-norm on both sides of (3.3), it holds that

Let $\phi(\mu_k) = \|R^{(k)}\|_{F}^2.$ The first-order derivative of $\phi(\mu_k)$ yields

It is easy to see that the unique stationary point of the function $\phi(\mu_k)$ is

It is obvious that the sencond-order derivative of $\phi(\mu_k)$ i.e., $\frac{\partial^2 \phi}{\partial \mu_k^2} = \frac{1}{2}\|M^{(k-1)}\|_{F}^2 > 0$, which implies that the stationary point mentioned in (3.4) is the unique minimum point of the function $\phi(\mu_k)$.

Through the above arrangement, we formally outline the APGI method in Algorithm 4.

Remark 5. If $P = I_{m}$ and $Q = I_{n}$ are adopted, the APGI algorithm is reduced to AGI algorithm.

3.2.   The adaptive GMI algorithm

The $k$th iteration of the GMI algorithm can be simply rewritten as:

The residual of the $k$$th$ iteration is

with $M^{(k-1)} = A^TR^{(k-1)}B+R^{(k-1)}B^TB+AA^TR^{(k-1)}+AR^{(k-1)}B^T,$ and $N^{(k-1)} = R^{(k-1)}-R^{(k-2)}.$ Taking the $F$-norm on both sides of (3.5), it follows that

Let $\psi(\mu_k, \beta_k) = \|R^{(k)}\|_{F}^2.$ The first-order derivative of $\psi(\mu_k, \beta_k)$ yields

It is easy to see that the unique stationary point of the function $\psi(\mu_k, \beta_k)$ is

where $a_{k-1} = tr((M^{(k-1)})^TR^{(k-1)})$, $b_{k-1} = tr((M^{(k-1)})^TN^{(k-1)})$, $c_{k-1} = tr((N^{(k-1)})^TR^{(k-1)})$, $d_{k-1} = \|M^{(k-1)}\|_F^2$, $e_{k-1} = \|N^{(k-1)}\|_F^2$. We also know that the sencond-order derivative of the function $\psi(\mu_k, \beta_k)$ is

whose Hessian matrix of the function $\psi(\mu_k, \beta_k)$ at the stationary point is

It is obvious that the Hessian matrix is symmetric positive definite, which implies that the stationary point mentioned in (3.6) is the unique minimum point of the function $\psi(\mu_k, \beta_k)$.

According to the above explanation, the AGMI algorithm is formulated as described in Algorithm 5.

4.   Numerical results

In this section, some examples are illustrated to verify the efficiencies of the proposed PGI, GMI, APGI, and AGMI algorithms compared with the GI ^[8], RGI ^[18], AGBI ^[25], AJGI ^[21], HSS ^[3], and NPHSS ^[16] algorithms. All examples are performed under MATLAB on a personal computer with a 1.61 GHz central processing unit (Intel(R) Core(TM) i7-10710) and 16GB memory.

The number of the iteration steps (denoted by IT), the computing time in seconds (denoted by CPU), and the relative residual norm (denoted by RRN) are listed in the tables below. All the initial matrices are set to be zero matrices, and the iterations are stopped if the RRN in the current step satisfies

or the numbers of iteration steps exceeds 10000.

Example 1. Consider Eq (1.1) with the matrices $A$ and $B$ defined by

with $L$ is the strictly lower triangular matrix having ones in the lower triangle part, $r = 2$ and $t = \frac{1}{2}$. The right-hand side $C = AX+XB$ with $X = ones(n)$, where $ones$ is a MATLAB built-in function.

The numerical results of the tested algorithms for Example 1 are listed in Table 2, and the corresponding error convergence curves are shown in Figure 1. The matrices $P$ and $Q$ in the PGI and APGI algorithms are the identity matrices, so the PGI and APGI algorithms are just the GI and AGI algorithms. The parameters in AGBI (Algorithm 2.4 in ^[25]), GI ((13)–(15) in ^[8]), RGI (Algorithm 1 in ^[18]), and GMI algorithms are experimentally optimal, which are denoted by $\mu_{exp}$ and $\beta_{exp}$ in Table 1. Like ^[18] and ^[25], the relaxation parameters in RGI and AJBI algorithms are both 0.5. Figure 1 shows the RRN of the AGBI, GI, AGI, GMI, and AGMI algorithms with $n = 200$. From the figure, it can be observed that the AGMI algorithm performs best among all these algorithms.

From Table 2, compared with the GI, RGI, AGBI, and AGI algorithms, the GMI and AGMI algorithms have more effectiveness in terms of IT and CPU time. In addition, since the parameters $\mu_k$ and $\beta_k$ in AGI and AGMI algorithms are varied and adaptive in each iteration, the two algorithms are more efficient than the GI and GMI algorithms, respectively.

Example 2. The matrices $A$ and $B$ in Eq (1.1) are given as

The right-hand side $C = AX+XB$ with $X = ones(n)$.

The numerical results of the tested algorithms for Example 2 are listed in Table 3, and the corresponding error convergence curves are shown in Figure 2. The matrices $P$ and $Q$ in the PGI and APGI algorithms are the diagonal parts of the matrices $A$ and $B$, respectively. The experimentally optimal parameters contained in GI (see (13)–(15) in ^[8]), HSS ("The HSS Iteration Method" in ^[3]), NPHSS (Algorithm 1 in ^[16]), PGI, and GMI algorithms are given in Table 4. From Table 3, it is clear that the seven algorithms are convergent for all cases. Furthermore, the APGI and AGMI algorithms need fewer iteration numbers and CPU time than other algorithms, and the AGMI algorithm performs best among the seven algorithms. Figure 2 illustrates the convergence curves of the HSS, NPHSS, PGI, GMI, APGI, and AGMI algorithms for the case $n = 256$. It is clear that the RRN of the AGMI and NPHSS algorithms rapidly decreases below $10^{-13}$, but NPHSS needs much more CPU time. So, AGMI is the most efficient one.

Example 3. The matrices in Eq (1.1) are described as

where $M = {\rm tridiag}(-1, 2.6, -1)$ and $N = {\rm tridiag}(0.5, 0, -0.5)$. The right-hand side $C = AX+XB$ with $X = ones(n)$.

The numerical results of the tested algorithms for Example 3 are listed in Table 6, and the corresponding error convergence curves are shown in Figure 3. Let $P$ and $Q$ be the tridiagonal parts of the matrices $A^{T}A$ and $BB^{T}$, respectively. The experimentally optimal parameters contained in GI ((13)–(15) in ^[8]), AJGI (Algorithm 5 in ^[21]), PGI, and GMI algorithms are given in Table 5. Like ^[21], the relaxation parameters $\omega_{1}$ and $\omega_{2}$ in the AJGI algorithm are 0.5 and 3, respectively. From Table 6 it follows that all the algorithms are effective for this example. In addition, the APGI algorithm is the most efficient one among the six algorithms. Figure 3 shows their convergence performances for the case $n = 256$, and the APGI algorithm has the remarkably best convergence result.

5.   Conclusions

In this paper, we provide two accelerated GI algorithms for solving Eq (1.1), which are the preconditioned gradient-based iteration (PGI) algorithm and the gradient-based momentum iteration (GMI) algorithm, respectively. Convergence analyses show that the proposed algorithms converge to the exact solution for any initial value with some assumptions. Moreover, the adaptive PGI and GMI algorithms are also established, and the adaptive parameters can be computed by minimizing the residual norms in the corresponding algorithms. Numerical experiments illustrate the excellent performances of our proposed algorithms. In addition, how to use the APGI and AGMI algorithms for solving other matrix equations will be investigated in our future work.

Use of Generative-AI tools declaration

The authors declare that they have not used Artificial Intelligence (AI) tools in the creation of this article.

Author contributions

Huiling Wang: gave the algorithms proposed in the manuscript, provided the numerical results and wrote the original draft of the manuscript; Nian-Ci Wu and Yufeng Nie: gave the clear guidance on the proof of the theorem and polished the language of the entire manuscript. All authors have read and agreed to the published version of the manuscript.

Acknowledgments

The work is supported by the National Natural Science Foundation of China (12201651), the Fundamental Research Funds for the Central Universities, South Central Minzu University (CZQ23004), and the Research Project Supported by Shanxi Scholarship Council of China(2023-117).

Conflict of interest

The authors declare no conflicts of interest.