Convergence analysis of a gradient iterative algorithm with optimal convergence factor for a generalized Sylvester-transpose matrix equation

1.   Introduction

It is well known that the fundamental of differential equations deals with the algebraic linear system

where $x(t)$ is an unknown vector-valued function and $A$ is a given square matrix. To analyse the stability of an equilibrium point of the system (1.1), it suffices to find a positive definite matrix $L$ such that $A^TL+LA$ is negative definite; see e.g., ^[1]. So, we need to solve the so-called Lyapunov equation

for some negative definite matrix $R$. To discuss the Eq (1.2), we investigate a more general form, namely, the generalized Sylvester-transpose matrix equation

where $A_i, B_i, C_j, D_j$ and $F$ are given matrices of conforming dimensions and $X$ is an unknown matrix to be determined. The Eq (1.3) includes important practical problems that are written as the matrix equations

known as Sylvester, Sylvester-transpose, generalized Sylvester, Kalman-Yakubovich matrix equations, respectively. The Eqs (1.2)–(1.7) have many essential applications in control theory; see e.g., ^[2,3,4,5,6]. In traditional method, a vectorization and the Kronecker product are used to find the unique solution. However, the large size of the Kronecker multiplication leads to computationally difficulty in that excessive computer. For this reason, iterative algorithms are received more attention.

Many researchers attempted to meliorate such iterative algorithms for finding the approximated solutions of matrix equations (1.2)–(1.7) using many ideas, e.g., matrix sign function ^[7,8], recursive blocked algorithms ^[9,10] and Hermitian and skew-Hermitian splitting algorithms ^[11,12,13]. In 2005, gradient-based iterative algorithms were firstly introduced by F. Ding and T. Chen for solving (1.4), (1.6) and (1.7); see ^[14,15]. In a few year later, many iterative algorithms that relied on gradients and hierarchical identification principle for solving (1.2)–(1.7) are established, e.g., RGI ^[16], JGI ^[17,18], MGI ^[19] and AGBI ^[20]. See more information in ^{[21,22,23,24,25]}. The Frobenious norm $\Vert\cdot\Vert_F$ and the spectral norm $\Vert\cdot\Vert_2$ for matrices are used to analyse the convergence property of such algorithms. There are defined respectively for any real matrix $A$ by

A gradient-based iterative algorithm for solving (1.3) was presented as follows:

Theorem 1.1 (^[25]). Suppose that the matrix equation (1.3) has a unique solution $X$. For each $s = 1, 2, \dots, p$ and $t = p+1, \dots, q$, construct

A necessary and sufficient condition for which the sequence $\{X_k\}$ converges to $X$ for any given initial matrices $X_1(0), X_2(0), \dots, X_{p+q}(0)$ is that

Meanwhile, a least-squares based iterative algorithm for solving (1.3) was presented as follows:

Theorem 1.2 (^[25]).

For each $s = 1, 2, \dots, p$ and $t = p+1, \dots, q$, construct

The sequence $\{X_k\}$ converges to a unique solution $X$ for any given initial matrices $X_1(0), X_2(0), \dots, X_{p+q}(0)$ if and only if $0 < \tau < 2(p+q)$.

In this work, we propose an effective iterative algorithm based on gradient and hierarchical identification principle, namely, a gradient iterative algorithm with optimal convergence factor. The proposed algorithm is applicable for the generalized Sylvester-transpose matrix equation (1.3); see Section 2. Convergence analysis (see Section 3) via Banach fixed-point theorem reveals that the iterative solutions converges to the unique solution for any initial value if and only if the convergence factor is chosen appropriately belong to an open interval. Then we study the performance of the algorithm from the convergence rate and error estimates. Moreover, we determine the fastest asymptotic convergence rate to minimize the spectral radius of the iteration matrix. Furthermore, we apply the algorithm to the Sylvester-transpose matrix equation; see Section 4. To show the applicability and the performance of the algorithm, we give numerical experiments in Section 5. In Section 6, we summarize the whole work. For benefits of reader, we include MATLAB-code for numerical experiment in Appendix.

2.   Introducing a gradient iterative algorithm

Let us denote by $\mathbb{R}^{r\times s}$ the set of $r\times s$ real matrices. Consider the matrix equation (1.3) where $A_i\in\mathbb{R}^{m\times n}$, $B_i\in \mathbb{R}^{p\times q}$, $C_j\in\mathbb{R}^{m\times p}$, $D_j\in\mathbb{R}^{n\times q}$, $F\in \mathbb{R}^{m\times q}$ are given coefficient matrices and $X\in \mathbb{R}^{n\times p}$ is an unknown matrix to be determined. The dimension matching of matrices is assumed to be $mq = np$.

2.1.   A traditional method for the generalize Sylvester-transpose matrix equation

Recall that the commutation matrix $K_{mn}$ is defined by

where $e_i^{n}$ is the $i$th column of the $n \times n$ identity matrix $I_n$. Essential properties of the commutation matrices are given in the following lemma:

Lemma 2.1. (see e.g., [26,Ch. 4]) For any $A\in\mathbb{R}^{m\times n}$ and $B \in \mathbb{R}^{p\times q}$, we have

A traditional method to solve the matrix equation (1.3) is to take the vector operator and utilize Lemma 2.1. Indeed, we arrive at an equivalent linear system $Qx\: = \:b$ where

$x = {\rm{vec}}[X]$ and $b = {\rm{vec}}[F]$. From now on, assume that $Q$ is invertible. So, the matrix equation (1.3) has a unique solution

However, if the dimensions of $A_i, B_i, C_j, D_j$ are not small, e.g., $10^2\times 10^2$, then the dimension of $Q$ is $10^4\times 10^4.$ Such a dimension problem leads to computational difficulty for computation and inversion of large matrices. Hence, this approach is only applicable for small dimensional matrices.

2.2.   Gradient iterative algorithm for the matrix equation

We now propose an effective iterative algorithm to solve (1.3).

If $p \geqslant q$ then we set $C_j = 0$ and $D_j = 0$ for any $j > q$. If $q > p$ then we set $A_i = 0$ and $B_i = 0$ for any $i > p$. For each $i = 1, 2, \ldots, p$, we assume that

From the main system (1.3), we would like to solve $p$ subsystems

so that the following errors are minimized:

We can derive the gradient of each $L_i$ as follows:

Let $X_i(k)$ be the approximated solution of the system (2.6) at iteration $k$. The recursive formula of $X_i(k)$ come from the gradient formula of $L_i(X)$ as follows:

where $\mu$ is a step size parameter. To avoid duplicated computation, we introduce a matrix

Taking the arithmetic mean of $X_1(k), \dots, X_p(k)$ to get $X(k)$:

where $\tau : = -2\mu/p$ is called a convergence factor. The hierarchical identification principle suggests to replace the unknown solution $X$ in (2.9) by its previous estimate $X(k-1)$. Thus we obtain the following procedure:

For convenience, we write $X^T(k)$ and $E^T(k)$ for $X(k)^T$ and $E(k)^T$, respectively. The matrices $E(k), A'_i, B'_i$ were introduced to avoid duplicate manipulations.

Remark 2.2. To break the procedure, if $\Vert F\Vert_F$ is close to zero, then we should consider the error $\Vert E(k)\Vert_F$ or $\lVert {X(k)-X(k-1)} \rVert_F$ instead of the relative error $\Vert E(k)\Vert_F/\Vert F\Vert_F$.

The convergent property of the algorithm relies on the convergence factor $\tau$, which will be discussed in the next section.

3.   Convergence analysis of GIO algorithm

In this section, we discuss a convergence criteria for the applicability of Algorithm 1. Then, we analyse the performance of the algorithm through its convergence rate and error estimates. The last job is to determine the convergence factor for which the algorithm fits with the fastest asymptotic behaviour. The main idea for the analysis starts with transforming a recursive equation of the error of approximated solutions into a fixed-point iteration $x(k) = Sx(k-1)$, where $S:\mathbb{R}^{np}\to \mathbb{R}^{np}$ is a contraction mapping. Then, we apply the following Banach fixed-point theorem to analysis the algorithm.

Theorem 3.1. (see e.g., ^[27]) Let $(X, d)$ be a non-empty complete metric space with a contraction mapping $T:X\to X$. Then

(i) The map $T$ admits a unique fixed-point $x^*$ in $X$ i.e. $T(x^*) = x^*$.

(ii) The fixed-point $x^*$ can be found as follows: start with an arbitrary element $x_0$ in $X$ and define a sequence $\lbrace x_n\rbrace$ by $x_n = T(x_{n-1})$ for $n\geq 1$. Then $x_n\to x^*$.

(iii) The following inequalities hold and describe the speed of convergence:

3.1.   Convergence criteria

From Algorithm 1, at each $k$-th iteration, we start with considering the error matrix $\hat{X}(k) = X(k)-X$. Indeed, we have

and

Thus, Lemma 2.1 implies that

Taking the vector operator to the Eq (3.4) and utilizing (3.5), we get

Letting $S = I_{np}-\tau Q^TQ$ and $x(k) = {\rm{vec}}\hat{X}(k-1)$, we get a linear iteration

Using Theorem 3.1 and [28,Theorem 1], we deduce that the following are equivalent:

(ⅰ) the sequence $\{X(k)\}$ converges to $X$ for any initial value $X(0)$;

(ⅱ) $S$ is a contraction mapping (here, we view $S:\mathbb{R}^{np}\to \mathbb{R}^{np}$ as a mapping);

(ⅲ) the spectral norm of $S$ is less than $1$.

Since $S$ is symmetric, all its eigenvalue are real. Note that the eigenvalues of $S$ are of the form $1-\tau\lambda$ where $\lambda$ is any eigenvalue of $Q^TQ$. We can compute the spectral radius of $S$ as follows:

Thus, $\rho[S] < 1$ if and only if

Since $Q$ is invertible, the matrix $Q^TQ$ is positive definite and hence, $\Vert Q\Vert_2^2 = \lambda_{\max}(Q^TQ) > 0$. The condition (3.8) now becomes

We summarize convergence criteria for Algorithm 1 as follows:

Theorem 3.2. Consider the matrix equation (1.3) under the assumption that the matrix $Q$ is invertible (i.e., Eq (1.3) has a unique solution). Let $\tau\in\mathbb{R}$. Then a necessary and sufficient condition for which Algorithm 1 is applicable for any initial matrix X(0) is the condition (3.9).

3.2.   Performance of the algorithm

We now discuss the performance of Algorithm 1 through its convergence rate and error estimates.

Considering Eq (3.6), we get

Since $S$ is symmetric, we have $\Vert S\Vert_2 = \rho[S]$. Thus the Eq (3.3) implies that

By induction, we obtain

According to the estimate (3.11), the asymptotic convergence rate of the algorithm relies on $\rho[S]$. Moreover, given an error $\epsilon > 0$, we would like to find the iteration number $k$ for which

By taking the 10th-base logarithms, we obtain the following equivalent requirement:

Moreover, by Theorem 3.1(ⅲ), we have

The following results are discussed to concluding the convergence rate and several estimates of the proposed algorithm.

Theorem 3.3. Suppose the convergence factor $\tau$ is chosen so that Algorithm 1 is applicable for any initial matrix X(0).

(i) The spectral radius $\rho[S]$ in (3.7) governs the asymptotic convergence rate of the algorithm.

(ii) Equations (3.10) and (3.11) show the error estimates $\Vert X(k)-X\Vert_F$ compared to the previous step and the first step, respectively. Meanwhile, the error at each iteration reduces from the previous one.

(iii) The prior and posterior estimates are presented in (3.14) and (3.15), respectively.

(iv) For each given error $\varepsilon > 0$, we have $\Vert X(k)-X\Vert_F < \epsilon$ after the $k$-th iteration for any $k\in\mathbb{N}$ that satisfies (3.13).

We can see that $A_i, B_i, C_j, D_j$ affect the convergence rate of Algorithm 1 but $E$ is not. However, the matrix $E$ is required for stopping process. Furthermore, If we take $\epsilon = 0.5\times 10^{-n}$ in (3.13) and $k$ satisfies

then the approximated $X(k)$ has an accuracy of $n$ decimal digit.

3.3.   Optimal convergence factor

The fastest convergence factor of Algorithm 1 is discussed intensively. Recall that the condition number of a matrix $A$ (relative to the spectral norm) is defined by

Assume that Eq (3.9) holds. The convergence rate of Algorithm 1 is the same as that of the linear iteration (3.6), and thus, it is given by the spectral radius (3.7) of the iteration matrix $S$. The fastest convergence rate is equivalent to the smallest of $\rho[S]$. Thus, we would like to minimize this spectral radius subject to the condition (3.9). Now, we apply the following lemma:

Lemma 3.4. (^[19]) For any real number $a, b$ with $b > a > 0,$ we have

The minimality is reached at $x_{opt} = \dfrac{2}{a+b}.$

Then the minimum value of $\rho[S]$ is

Meanwhile, we can notice that spectral radius of the iteration matrix is

Thus, we obtain the following theorem:

Theorem 3.5. Among the convergence factors $\tau$ that meet the criteria of Algorithm 1, the one in Eq (3.16) attains the fastest asymptotic convergence rate of the algorithm, which is governed by the spectral radius (3.17).

The above theorem tells us that if $\lambda_{\max}(Q^TQ)$ is neighbouring to $\lambda_{\min}(Q^TQ)$, or equivalently, the condition number is close to one then Algorithm 1 has a fast convergence.

4.   The proposed algorithm for Sylvester-transpose matrix equation

In this section, we consider an important special case of the matrix equation (1.3), namely, the Sylvester-transpose matrix equation. We apply GIO algorithm for this equation, and investigate the convergence property of the algorithm.

Let $m, n\in\mathbb{N}$. Consider the Sylvester-transpose matrix equation (1.5) where $A\in\mathbb{R}^{m\times n}$, $B\in\mathbb{R}^{n\times m}$, $F\in \mathbb{R}^{m\times m}$ are given constant matrices and $X\in \mathbb{R}^{n\times m}$ is an unknown matrix to be solved. Suppose that (1.5) has a unique solution, i.e., that following matrix is invertible:

Corollary 4.1. Consider the matrix equation (1.5) when the matrix $P$ is invertible. Let $\tau\in\mathbb{R}$. Then the given statements hold:

(i) A equivalent condition for which Algorithm 2 is applicable for any initial matrix $X(0)$ is

Meanwhile, the spectral radius of the iteration matrix $T = I_{np}-\tau P^TP$ is given by

(ii) The spectral radius $\rho[T]$ in (4.2) represents the asymptotic convergence rate of Algorithm 2.

(iii) The fastest asymptotic convergence factor is determined by

5.   Numerical simulations with discussion

In this section, we report some numerical results to illustrate the applicability the effectiveness of Algorithm 1. All simulations have been carried out by MATLAB R2018a, AMD Ryzen7 3700U with Redeon Vega Mobile Gfx @ 2.30 GHz, RAM 12.00 GB PC environment. In Example 5.1, we show that our algorithm is applicable for small-square-matrices in the case $p > q$. We also show that our algorithm is applicable and efficient for large-square-matrices and consider the effect of changing the convergence factor $\tau$ in Example 5.2. In Example 5.3, we illustrate the efficiency of Algorithm 1 when coefficients are non-square matrices of different moderate sizes. In Examples 5.4 and 5.5, we do experiments in the cases $p > q$ and $q > p$ with rectangular coefficient matrices. In Example 5.6, we test Algorithm 2 for the Sylvester-transpose equation with square coefficient matrices. In all examples, we compare the proposed algorithm to both the traditional method (Eq (2.4)) and recent iterative algorithms. The computational time is measured in seconds by MATLAB functions tic and toc.

Example 5.1. Consider the matrix equation

where

In fact, the unique solution is given by

Let us apply Algorithm 1 to compute the sequence $\lbrace X(k)\rbrace$ of approximated solutions. Take an initial point $X(0) = 0.$ The optimal convergence factor can be computed according to Theorem 3.2 as follows:

Table 1 shows that the direct method consumes 15 ms to get the exact solution, while GIO algorithm takes 6 ms to perform $10$ iterations in order to get an approximated solution with a small relative error. Figure 1 and Table 1 illustrate that the computational time of GIO algorithm is less than that for GI algorithm with preferable relative error. Figure 1 and Table 1 also show that GIO algorithm is outperform than LS algorithm.

Example 5.2. Consider the matrix equation

where all coefficients are $100\times 100$ tridiagonal matrices given by

In fact, the unique solution is given by $X = {\rm{tridiag}}(0, 1, -1)$. Let us apply Algorithm 1 to compute the sequence $\lbrace X(k)\rbrace$ of approximated solutions. Take an initial point

The optimal convergence factor can be computed according to Theorem 3.2 as follows:

The computing results of changing $\tau$ are listed in Figure 2 and Table 2. Figure 2 shows that, for the given initial point, as $k$ increases, for $\tau = \tau_{opt}, \tau = 0.00002, \tau = 0.000015,$ and $\tau = 0.00001,$ we see that the terms $\lVert {E(k)} \rVert_F / \lVert {F} \rVert_F$ are becoming smaller and approaching to zero. Moreover, the relative errors $\lVert {E(k)} \rVert_F / \lVert {F} \rVert_F$ for $\tau_{opt}$ go faster to $0$ than those for another convergence factors. Furthermore, when $\tau = 0.00004$ and $\tau = -0.00001$ which do not satisfy the criteria (3.9), the approximated solutions diverge. Table 2 confirms that the computational time of the proposed algorithm is significantly less than the time of the traditional method. Thus, in this case, Algorithm 1 is applicable and effective.

Example 5.3. Consider the matrix equation

when $A_i \in \mathbb{R}^{40\times 60}$, $B_i \in \mathbb{R}^{20\times 30}$, $C_j \in \mathbb{R}^{40\times 20}$, $D_j \in \mathbb{R}^{60\times 30}$ and $F\in\mathbb{R}^{40\times 30}$ are tridiagonal matrices given by

Then the unique solution is $X = {\rm{tridiag}}(0, 1, -1) \in \mathbb{R}^{60 \times 20}$. We take the initial matrix

The computing results of changing $\tau$ are listed in Figure 3. As $k$ large enough, the terms $\lVert {E(k)} \rVert_F / \lVert {F} \rVert_F$ for $\tau_{opt} \approx 0.000085$ are becoming smaller and go faster to zero than another convergence factors. Moreover, for $\tau = 0.0003$ and $\tau = -0.000001$ which do not satisfy the condition (3.9), we see that the term $\lVert {E(k)} \rVert_F / \lVert {F} \rVert_F$ do not approach zero, so the approximated solutions diverge. From Table 3, it is clear that the computational time of GIO algorithm is significantly less than the time of the traditional method and another GI with different convergence factors.

Example 5.4. Consider the matrix equation

when $A_i \in \mathbb{R}^{40\times 60}$, $B_i \in \mathbb{R}^{20\times 30}$, $C_j \in \mathbb{R}^{40\times 20}$, $D_j \in \mathbb{R}^{60\times 30}$ and $F\in\mathbb{R}^{40\times 30}$ are tridiagonal matrices given by

Then the unique solution is $X = {\rm{tridiag}}(1, 2, -1) \in \mathbb{R}^{60 \times 20}$. We take the initial matrix

The numerical results in Figure 4 show that the direct method consumes around $11$ seconds, while GIO algorithm spends $0.047$ seconds ($50$ iterations) with satisfactory error. Figure 4 and Table 4 show that the computational time of GIO algorithm is slightly more than that of GI algorithm, as the relative error of GIO is significantly less than that of GI algorithm. Figure 4 and Table 4 also show that GIO algorithm is outperform than LS algorithm in both computational time and relative error.

Example 5.5. Consider the matrix equation

where all coefficients are $10\times 10$ tridiagonal matrices given by

In fact, the unique solution is given by $X = {\rm{tridiag}}(1, 1, 1)$. Let us apply Algorithm 1 to compute the sequence $\lbrace X(k)\rbrace$ of approximated solutions using an initial point

Table 5 shows that for small size matrices, the difference between the computational times for the direct method and GIO algorithm ($50$ iterations) is not much as that for the moderate sizes in Example 5.4. Figure 5 and Table 5 illustrate that the computational time of GIO algorithm is less than that for GI algorithm with preferable relative error. Figure 5 and Table 5 also show that GIO algorithm is outperform than LS algorithm.

Example 5.6. The Sylvester-transpose equation $AX+X^TB = F$, where $A, B, F, X\in \mathbb{R}^{10\times 10}$ are given by

is considered by using the initial matrix $X(0) = {\rm{tridiag}}(10^{-6}, 10^{-6}, 10^{-6})$. The objective is to compare the capability of GIO algorithm to GI ^[25], LS ^[25], and AGBI ^[20] algorithms. We fix the iteration number to be $50$ and investigate the relative error $\Vert E(k)\Vert_F/ \Vert F\Vert_F$. Figure 6 and Table 6 indicate that Algorithm 2 is more efficient than the traditional method and another iterative algorithms in spite of a little more computational time.

6.   Conclusions

The proposed algorithm, a gradient iterative algorithm with an optimal convergence factor, is applicable for a generalized Sylvester-transpose matrix equation (1.3) under the assumption that the matrix equation has a unique solution. The analysis tells us that a necessary and sufficient condition for which Algorithms 1 and 2 are applicable for any initial matrix $X(0)$ is the convergence factor is chosen appropriately belong to an open interval. The convergence rate and several estimations (e.g., error, prior, posterior) depend on the spectral radius of the associated iteration matrix. Moreover, we can determine the fastest asymptotic convergence rate. The numerical experiments illustrate that the algorithm can be applied for any matrix problems with conformable coefficient matrices of small/moderate/large sizes and square/non-square sizes. Moreover, the algorithm has more efficient than the traditional method and another recent gradient-based iterative algorithms.

Acknowledgements

The first author would like to thank Ministry of Education, Thailand for a financial support during her PhD study.

Conflict of interest

The authors declare that they have no conflict of interest.

Appendix

This appendix contains MATLAB-code of Example 5.1

tic;

$\tau = 0.1379;$

Xt = X;

XtT = transpose(Xt);

Et = F - (A1*Xt*B1 + A2*Xt*B2 + A3*Xt*B3 + C1*XtT*D1 + C2*XtT*D2 + C3*XT*D3);

EtT = transpose(Et);

et(1) = e(1)/norm(F, 'fro');

let(1) = le(1);

kt = 2;

while $(kt < 50)$

Et = F - (A1*Xt*B1 + A2*Xt*B2 + A3*Xt*B3 + C1*XtT*D1 + C2*XtT*D2 + C3*XT*D3);

EtT = transpose(Et); Xt = Xt + $\tau$*(A1T*Et*B1T+ A2T*Et*B2T + A3T*Et*B3T + D1*EtT*C1 + D2*EtT*C2 + D3*EtT*C3);

et(kt) = norm(Et, 'fro')/norm(F, 'fro');

let(kt) = log(norm(Et, 'fro'));

kt = kt+1;

end

tt = toc;