Parametric Smith iterative algorithms for discrete Lyapunov matrix equations

An iterative algorithm is established in this paper for solving the discrete Lyapunov matrix equations. The proposed algorithm contains a tunable parameter, and includes the Smith iteration as a special case, and thus is called the parametric Smith iterative algorithm. Some convergence conditions are developed for the proposed parametric Smith iterative algorithm. Moreover, the optimal parameter for the proposed algorithm to have the fastest convergence rate is also provided for a special case. Finally, numerical examples are employed to illustrate the effectiveness of the proposed algorithm.


1.
Introduction. Discrete Lyapunov matrix equations in the form of AXA T −X = −Q play important roles in analysis and design of discrete-time linear systems [7]. A well-known result is that the stability of a discrete-time linear system can be characterized by the existence of the corresponding discrete Lyapunov matrix equation. In model reduction of discrete-time systems, impulse response Gramians are very important, and they can be used to derive the reduced-order model. In [13] and [1], the concepts of weighted and generalized impulse response Gramians were respectively proposed, and it was shown that these Gramians can be obtained by solving the corresponding discrete Lyapunov matrix equations. In addition, the discrete Lyapunov matrix equation can also be used to analyze the stability of periodic discrete-time systems [12].
Due to the wide application of discrete Lyapunov matrix equations, many researchers have made much effort on the solution of this kind of matrix equations. In [11], a discrete Lyapunov matrix equation in the controllable canonical form was investigated. It was shown that the solution is inverse of the corresponding Schur-Cohn matrix. In [6], an algorithm for obtaining the solution of the discrete Lyapunov matrix equation AXA T − X = −Q was provided in a factored form. In this algorithm, real Schur decomposition for the matrix A is needed, and the original matrix equation is transformed into a quasi-triangular Lyapunov matrix equation. In [17], a new updating formula was established for the Cholesky factors in the Hammarling's algorithm of [6]. In [5], an improved rank 2 updating formula was proposed for the Cholesky factors arising in the numerical algorithm to solve the discrete Lyapunov matrix equation. In [15], an iterative algorithm was given to solve the discrete Lyapunov matrix equation by means of Kronecker products. Obviously, a higher-dimensional matrix is involved when the algorihtm in [15] is used. A common feature of the algorithms in [5], [6], [11], [15], [17] is that the coefficient matrices of the equations need to be transformed into or be in some canonical forms.
In addition, there have existed some algorithms which are implemented by only using original coefficient matrices. For Kalman-Yakubovich matrix equation AXB− X = C, there exists the well-known Smith iteration in the form of X(m + 1) = AX(m)B + C [10]. Obviously, this algorithm can also be applied to solve the discrete Lyapunov matrix equation. In [4], a gradient based iterative algorithm was proposed to solve the matrix equation AXB − X = C by applying the hierarchical identification principle, and the convergence condition was also given by original coefficient matrices. By specializing the algorithm in [4], a gradient based iterative algorithm can also be obtained to solve the discrete Lyapunov matrix equation. It should be pointed out that the approaches in [10] and [4] are different from those in [5], [6], [11], [15], and [17]. In [10] and [4], it is not necessary to transform the coefficient matrices into any canonical form.
In addition, it should be pointed out that the dual version AX − XB = C of the matrix equation AXB − X = C was also extensively investigated. In [4], by using the hierarchical principle a gradient based iterative algorithm was proposed. A modified version of the algorithm in [4] was presented in [9] to solve the matrix equation AX − XB = C. In this algorithm, each iterative step needs to be divided into two substeps, and thus a so-called relaxation parameter was introduced. In [16], the matrices A and B were firstly written as the sum of diagonal matrices and other matrices, and then a Jacobi-like algorithm was proposed. Recently, the Lyapunov matrix equations appearing in stochastic linear systems were investigated in [18] and [19], and iterative algorithms were presented to solve this kind of matrix equations in the case where the associated systems are stochastically stable. In addition, a kind of coupled matrix equations in a general form was investigated in [3], and an iterative method was developed to solve this kind of matrices by extending the idea of conjugate gradient methods.
In this paper, we are interested in the iterative algorithms constructed by using the original matrices of the discrete Lyapunov matrix equations. By introducing a tuning parameter, an equivalent expression of the considered discrete Lyapunov matrix equation is first given, and then an iterative algorithm is constructed. The proposed algorithm includes the Smith iteration as a special case, and thus is called parametric Smith iteration in this paper. Some convergence conditions are given for the parametric Smith iterative algorithm. The motivation of this paper lies two aspects. one is that the original Smith iterative algorithm is only valid for the case that the corresponding system is Schur stable. In this paper, it will be seen that the algorithm in the current paper can be applicable for some discrete Lyapunov matrix equations with unstable matrices.
Throughout this paper, we use · 2 , · F and ρ(·) to denote the spectral norm, Frobenius norm and spectral radius of a matrix, respectively. In addition, σ(·) is used to denote the eigenvalue set of a matrix. For two integers a ≤ b, the notation I[a, b] is defined as I[a, b] = {a, a + 1, ..., b}. For a complex number µ, we useμ to denote its conjugate.

2.
Previous results and parametric Smith iterative algorithms. In this paper, we consider the following discrete Lyapunov matrix equation where A ∈ R n×n and Q ≥ 0 are known matrices, and X ∈ R n×n is the matrix to be determined. This matrix equation plays a vital role in analysis and design of discrete-time linear systems. For example, it can be used to check the stability of a discrete-time linear system. In addition, the solution of this kind of matrix equations can be used to characterize the H 2 performance, controllability and observability. As in [10], with the special structure of the discrete Lyapunov matrix equation (1) the following Smith iterative algorithm is easily obtained: In addition, by specializing the gradient based iterative algorithm in [4] the following result can be obtained to solve the discrete Lyapunov matrix equation (1).
converges to the unique solution of (1) if the parameter µ satisfies In this section, we aim to give another iterative algorithm for solving the discrete Lyapunov matrix equation (1). For the unknown matrix X in the discrete Lyapunov matrix equation (1), for any γ ∈ R there holds Combining this expression with (1), gives Based on this identity, the following iterative algorithm can be constructed to solve the discrete Lyapunov matrix equation (1): If γ = 1, then the algorithm (5) becomes the Smith iteration in (2). Due to this reason, for convenience the algorithm (5) is called parametric Smith iteration.
3. Monotonicity and boundedness. In this section, monotonicity and boundedness are considered for the proposed parametric Smith iterative algorithm (5). First, let us show monotonicity of the sequence generated by the algorithm (5).
Lemma 3.1. Given a matrix A to be Schur stable, the sequence {X (m)} generated by the algorithm (5) with γ ∈ (0, 1] for given positive definite matrices Q is strictly monotonically increasing under zero initial condition X(0) = 0. That is, for any integer m ≥ 0, there holds X (m + 1) > X (m) .
Proof. We use the principle of mathematical induction to prove the conclusion of this lemma. First, due to the zero initial condition, it is obvious that (6) is true for m = 0. Now, it is assumed that (6) holds for m = s, s ≥ 0, that is With the induction assumption (7), it follows from the expression in (5) with m = s + 1 and γ ∈ (0, 1] that This implies that (6) holds for m = s + 1. The previous relations (8) and (6) with m = 0 imply that (6) holds for any m ≥ 0 by mathematical induction. The proof is thus completed.
Next, the boundedness of the sequence generated by the algorithm (5) is investigated. The result is given in the following lemma.
Lemma 3.2. Given a Schur stable A, the sequence {X (m)} generated by the algorithm (5) with γ ∈ (0, 1] for a given positive definite matrix Q is upper bounded by the solution of the discrete Lyapunov matrix equation (1) under zero initial condition X(0) = 0. That is, for any integer m ≥ 0, there holds Proof. Similarly, we use the principle of mathematical induction to prove the conclusion of this lemma. Due to the zero initial condition, the expression (9) holds for m = 0. Now it is assumed that (9) holds for m = s, s ≥ 0, that is With the induction assumption (10), it follows from the expressions (4) and (5) with m = s and γ ∈ (0, 1] that This implies that the inequality (9) is true for m = s + 1. The relations (11) and (9) with m = 0 imply that (9) holds for any m ≥ 0 by mathematical induction principle.
With the monotonicity and boundedness properties given in the previous lemmas, a convergence result of the algorithm (5) is given in the following theorem. Proof. According to Lemmas 3.1 and 3.2, the sequence {X(m)} generated by the parametric Smith iterative algorithm (5) is monotonically increasing and upper bounded for γ ∈ (0, 1]. Thus, by the convergence principle of positive definite matrices in [2] the sequence {X(m)} is convergent. Denote lim m→∞ X(m) = X * . Then, by taking limits for (5), we have Since γ = 0, it follows from the preceding relation that This implies that X * is a solution of the discrete Lyapunov matrix equation (1).
Since the matrix A is Schur stable, then the Lyapunov matrix equation (1) has a unique positive definite solution. Thus, X * is the unique solution of the matrix equation (1). The proof is thus completed.
4. Convergence conditions. In the preceding section, the properties of monotonicity and boundedness of the sequence generated by the algorithm (5) are investigated. In this case, the following three requirements are needed.
In the section, we further investigate the convergence property of the algorithm (5) without the requirement of Schur stability and the zero initial condition. Moreover, the tuning parameter γ is not required to be confined to the interval (0, 1]. Theorem 4.1. Consider the discrete Lyapunov matrix equation (1) to have a unique solution, and define the following matrix The sequence {X (m)} generated by the algorithm (5) under an arbitrary initial condition converges to the solution of the discrete Lyapunov matrix equation (1) if and only if the parameter γ satisfies ρ(G) < 1.
Proof. Taking the vectorization for both sides of (5), gives vec (X (m + 1)) (14) Then, the relation (14) can be written as By the stability theory of discrete-time linear systems, the iteration (15) is convergent if and only if the parameter γ satisfies ρ(G) < 1. Further, let η = vec (X), then the discrete Lyapunov matrix equation (1) can be equivalently expressed as (I − A ⊗ A)η = q, and thus its unique solution can be given by η = (I − A ⊗ A) −1 q. In addition, it is easily derived that the vector η (m) in (15) converges to γ(I − G) −1 q under the condition of ρ(G) < 1. By using the expression in (13) one has This implies that the conclusion of this theorem is true.
In Theorem 4.1, the convergence condition for the algorithm (5) is given in terms of the spectral radius of a matrix with the tuning parameter γ. In the following theorem, an explicit necessary conditon is provided for the convergence of the algorithm (5).
Theorem 4.2. Let A be given, and µ i , i ∈ I[1, n], be the n eigenvalues of the matrix A, and denote The sequence {X (m)} given by the algorithm (5) converges to the unique solution of the discrete Lyapunov matrix equation (1) if and only if the parameter γ is chosen such that for any i, j ∈ I[1, n], there holds |λ ij (γ)| < 1.
Proof. Let λ be an arbitrary eigenvalue of the matrix G in (13), and v be the corresponding eigenvector. Thus, we have Gv = λv.
With the expression of (13), it follows from this relation that It is easily known that γ = 0 if the algorithm (5) is convergent. Thus, from (17) we have This implies that λ+γ−1 γ is an eigenvalue of A ⊗ A if λ is an eigenvalue of G. In addition, by properties of Kronecker products it is derived that From this, it can be derived that all the n 2 eigenvalues λ ij (γ) of the matrix G for a fixed γ are given by (16). Combining this with Theorem 4.1, gives the conclusion of this theorem.
By using Theorem 4.2, a convergence condition is established for the algorithm (5) in the following theorem.
Then, there exists a parameter γ such that the parametric Smith iterative algorithm (5) converges to the unique solution of the discrete Lyapunov matrix equation (1) if and only if if the condition (20) is satisfied, then the parametric Smith iterative algorithm (5) converges to the unique solution of the discrete Lyapunov matrix equation (1) if and only if Proof. According to Theorem 4.2, the algorithm (5) converges to the unique solution of the discrete Lyapunov matrix equation (1) if and only if for any i ∈ [1, s], the parameter γ satisfies |1 − γ + γ(a i ± ib i )| < 1, which is equivalent to (1 − γ + γa i ) 2 + (γb i ) 2 < 1. By simple calculation, from the preceding expression we have (a i − 1) 2 + b 2 i γ + 2(a i − 1) γ < 0. From this relation, the conclusion of this theorem can be immediately obtained.
Remark 1. In Theorem 3.3, the convergence condition of the algorithm (5) is given under the requirement of the Schur stability of the matrix A. However, the convergence conditions given in Theorems 4.1-4.3 are for a general discrete Lyapunov matrix equation in the form of (1) without the requirement that the matrix A is Schur stable. (2) is only valid for the discrete Lyapunov matrix equation (1) with A Schur stable. However, the proposed parametric Smith iteration can be used to solve a discrete Lyapunov matrix equation without the Schur stability condition for A if the tuning parameter is properly chosen.

Remark 2. It has been known that the Smith iteration
For the discrete Lyapunov matrix equation (1), denote the eigenvalue set of A⊗A as in (18). If the matrix A is Schur stable, then |a i ± ib i | < 1 for any i ∈ I[1, s]. Thus, |a i | < 1 for any i ∈ I [1, s]. With this, by Theorem 4.3 we can obtain the following result. In addition, another proof is also given for this result by using Theorem 4.2. Proof. Let µ be an eigenvalue of A. Then µ is also an eigenvalue of A. Since A is Schur stable, then 0 ≤ µµ = |µ| 2 < 1. Now, it is assumed that γ ≤ 0, then According to Theorem 4.2, the algorithm (5) is not convergent. This contradicts the convergence of the algorithm (5). Therefore, the assumption γ ≤ 0 is not true. Thus, there holds γ > 0.
A simple necessary condition is given in Lemma 4.4 for the algorithm (5) to be convergent when the matrix A is Schur stable. Such a result will play a vital role in the next section.

5.
Optimal tuning parameters. In the previous section, some convergence conditions have been established for the parametric Smith iterative algorithm (5). In this section, we focus on the choice of the tuning parameter such that the algorithm (5) has the fastest convergence rate. Firstly, we consider the case where the matrix A is Schur stable. {µ i µ j } < 1, τ = min i,j∈I [1,n] {µ i µ j } > −1, then the algorithm (5) under an arbitrary initial condition has the fastest convergence rate when the parameter is chosen as In this case, the spectral radius ρ(G) of the matrix G given in (13) is as follows Proof. For the algorithm (5), define the matrix G as in (13). From Theorem 4.2, the eigenvalue set of the matrix G is According to Theorem 4.4, under the condition of this theorem the parameter γ should be chosen to satisfy γ < 0 if the algorithm (5) is convergent. For some µ i µ j ≥ γ−1 γ , there holds 1 − γ + γµ i µ j ≥ 0. In this case, we have For some µ i µ j < γ−1 γ , there holds 1 − γ + γµ i µ j < 0. In this case, we have Now, we aim to give the expression of the spectral radius of the matrix G by using the previous preliminary. If κ < γ−1 γ , that is, γ > 1 1−κ , then we have If the parameter γ satisfies τ < γ−1 γ ≤ κ, that is, Now, we consider the case where τ ≥ γ−1 γ . In this case, the parameter γ satisfies With the previous derivation, it can be obtained that It is easily derived that By taking the monotonicity of the function in (26) into consideration, it is known that the spectral radius ρ(G) achieves its minimum when the parameter γ is chosen as in (21). In this case, the spectral radius ρ(G) is as in (22).

Remark 3.
In the proof of Theorem 5.1, the expression (26) of the spectral radius of the matrix G defined in (13) has been provided when the matrix A is Schur stable and all the eigenvalues of A are real. In view of the monotonicity of the function in (26), it is easily known that in this case the algorithm (5) is convergent if and only if the following conditions are satisfied which is equivalent to This condition can also be obtained from Theorem 4.3. In fact, if all the eigenvalues of A are real, then b i = 0, a i < 1, i ∈ I [1, s]. In this case, the convergence condition in Theorem 4.3 is reduced to This condition is equivalent to the condition (27).
In the sequel, we consider the algorithm (5) for the discrete Lyapunov matrix equation (1) where all the eigenvalues of A are real and greater than 1. According to Theorem 4.3, a necessary condition for the algorithm (5) to be convergent is γ < 0.
Theorem 5.2. Consider the discrete Lyapunov matrix equation (1) to have a unique solution, and let µ i , i ∈ I[1, n], be the n eigenvalues of the matrix A. If all µ i µ j , i, j ∈ I[1, n], are real and τ = min i,j∈I [1,n] {µ i µ j } > 1, then the algorithm (5) under an arbitrary initial condition has the fastest convergence rate when the parameter is chosen as where κ is defined as in Theorem 5.1. In this case, the spectral radius ρ(G) of the matrix G given in (13) is as follows Proof. From Theorem 4.2, the eigenvalue set of the matrix G is According to Lemma 4.4, under the condition of this theorem the parameter γ should be chosen to satisfy γ > 0 if the algorithm (5) is convergent. For some µ i µ j ≤ γ−1 γ , there holds 1 − γ + γµ i µ j ≥ 0. In this case, we have For some µ i µ j > γ−1 γ , there holds 1 − γ + γµ i µ j < 0. In this case, we have With the preceding preliminary, in the sequel we give the expression of the spectral radius of the matrix G. If κ < γ−1 γ , that is, γ < 1 1−κ , then we have ρ(G) = max i,j∈I [1,n] {|1 − γ + γµ i µ j |} = max i,j∈I [1,n] {1 − γ + γµ i µ j } If the parameter γ satisfies τ < γ−1 γ ≤ κ, that is, 1 1−κ ≤ γ < 1 1−τ , then we have Now, we consider the case where τ ≥ γ−1 γ . In this case, the parameter γ satisfies γ ≥ 1 1−τ . We have With the previous derivation, it can be obtained that It is easily derived that With this, it follows from (32) that The conclusion of this theorem can be obtained by using the similar treatment in Theorem 5.1.

Remark 4.
In Theorems 5 and 6, not only the optimal parameter of the proposed parametric Smith iterative algorithm (5) is given, but also the corresponding spectral radius of the convergence matrix is provided. In addition, it can be seen from the proof that the spectral radius of the convergence matrix is a linear function with respect to the tuning parameter.
Remark 5. In this paper, the optimal tuning parameter of the proposed parametric Smith iterative algorithm (5) is provided only for the matrix A to have real eigenvalues. For the case where the matrix A has complex eigenvalues, it is difficult to obtain an explicit expression of the optimal tuning parameter.
and Q = I. For this example, the Smith iterative algorithm (2) does not work. By trial and test, it is found that the gradient based algorithm (3) achieves its fastest convergence rate when µ = 8.12 × 10 −3 . Denote The convergence curve of the gradient based algorithm (3) for this example is given in Fig. 1 when the initial value is chosen to be X(0) = 0. When the proposed algorithm (5) is applied to this example, it is found by trial and test that the sequence generated by this algorithm is convergent when −0.13 < γ < 0. In fact, this convergence range can be obtained by using Theorem 4.3. In addition, the spectral radius of G defined in (13) versus the tuning parameter γ is shown in Fig. 2. According to this figure, the spectral radius of G is minimum when the parameter γ is chosen as γ = −0.1279. Therefore, the proposed algorithm (5) achieves its fastest convergence rate when γ = −0.1279. Shown in Fig. 3 is the convergence curve of the proposed algorithm (5) for different parameters γ under the zero initial condition, that is, X(0) = 0. It follows from Figures 1 and 3 that the proposed parametric Smith iterative algorithm (5) converges much faster than the gradient algorithm (3).   (5) is shown in Fig. 6 for this example. It should be pointed out the parametric Smith iterative algorithm (5) with γ = 1 is the Smith iteration algorithm. It can be seen from Figures 4 and 6 that the proposed parametric Smith iterative algorithm has better convergence performance than the existing Smith iterative algorithm and gradient based algorithm if the tuning parameter is appropriately chosen. Thus, for this example κ = 0.5625 and τ = −0.3525. According to the result of Theorem 5.1, the parametric Smith iterative algorithm (5) with 0 < γ < 1.478 is convergent, and the optimal tuning parameter is γ = 1.117. In addition, the spectral radius of G for this example is given in Fig. 7. It is easily found that the spectral radius of G achieves its minimum when γ = 1.117. Moreover, the convergence curves for different tuning parameters are given in Fig. 8 when the initial matrix X(0) is chosen as From this figure, it can be seen that the algorithm (5) has the fastest convergence rate when γ = 1.117. These facts illustrate the effectiveness of Theorem 5.1.

7.
Conclusions. In this paper, we have developed a novel iterative algorithm to solve the well-known discrete Lyapunov matrix equations. This algorithm is established by introducing a tunable parameter, and includes the well-known Smith iterative algorithm as a special case. The convergence properties of this kind of algorithms are analyzed, and some convergence conditions are provided. In addition, a choice approach is also established for the optimal tuning parameter of the proposed algorithm. It has been shown that the proposed parametric Smith iterative algorithm could have better convergence rate than some existing iterative algorithms if the tuning parameter is properly chosen. Recently, a finite iterative algorithm was proposed in [8] to solve a class periodic Sylvester matrix equations. In [14], an inner-outer iterative algorithm was presented to solve the coupled continuous Markovian jump Lyapunov matrix equations. It should be pointed out that the idea in the current paper can be used to solve these two kinds of matrix equations.