- Email: [email protected]

A property of the eigenvalues of the symmetric positive definite matrix and the iterative algorithm for coupled sylvester matrix equations Huamin Zhang, Feng Ding

www.elsevier.com/locate/jfranklin

PII: DOI: Reference:

S0016-0032(13)00313-X http://dx.doi.org/10.1016/j.jfranklin.2013.08.023 FI1860

To appear in:

Journal of the Franklin Institute

Received date: Revised date: Accepted date:

4 February 2013 23 August 2013 25 August 2013

Cite this article as: Huamin Zhang, Feng Ding, A property of the eigenvalues of the symmetric positive definite matrix and the iterative algorithm for coupled sylvester matrix equations, Journal of the Franklin Institute, http://dx.doi.org/10.1016/j.jfranklin.2013.08.023 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

A property of the eigenvalues of the symmetric positive definite matrix and the iterative algorithm for coupled Sylvester matrix equations Huamin Zhanga,b , Feng Ding∗,a a Key

Laboratory of Advanced Process Control for Light Industry (Ministry of Education), Jiangnan University, Wuxi 214122, PR China b Department of Mathematics and Physics, Bengbu College, Bengbu 233030, PR China

Abstract In this paper, we discuss the properties of the eigenvalues related to the symmetric positive deﬁnite matrices. Several new results are established to express the structures and bounds of the eigenvalues. Using these results, a family of iterative algorithms are presented for the matrix equation AX=F and the coupled Sylvester matrix equations. The analysis shows that the iterative solutions given by the least squares based iterative algorithms converge to their true values for any initial conditions. The eﬀectiveness of the proposed iterative algorithm is illustrated by a numerical example. Key words: Eigenvalue, Iterative algorithm, Matrix equation, Symmetric positive deﬁnite matrix

1. Introduction Matrix equations are often encountered in the system theory, control theory and some areas of the pure and applied mathematics [1, 2, 3]. Many publications have studied how to solve diﬀerent types of matrix equations, e.g., [4]. Though the formula solutions are vital in the theoretical derivation, the numerical solutions play the important roles in the practical applications [5, 6, 7]. If the dimensions of the coeﬃcient matrices are small then the Gaussian elimination or other direct methods are eﬀective. With the increase of the sizes of the related matrices, the iterative approaches have replaced the direct methods and become the main strategy of solving the matrix equations [8]. So the iterative algorithms for ﬁnding the numerical solutions have received mach attentions in many areas [9, 10], including signal processing and parameter estimation [11, 12, 13, 14] and system modeling [15, 16, 17, 18, 19]. Recently, by using the hierarchical identiﬁcation principle, Ding and his coworkers presented several iterative schemes for solving matrix equations [20, 21, 22]. The schemes, resembling the classical Jacobi and Gaussian iterations [5, 6] for linear systems, are easy to implement and cost little per step and are convergent linearly at the best. Gradient and least squares based iterative algorithms have also been proposed for (coupled) matrix equations [23, 24], which are applicable to the Lyapunov matrix equations and Sylvester matrix equations as special cases [25]. Inspired by this work, with the real representations of complex matrices, Wu et al. presented the gradient based iterative algorithms for solving the complex (coupled) matrix equations [26]; Song et al. proposed the gradient based iterative algorithms for the generalized coupled conjugate transpose matrix equations [27]. Contrasting to the ﬂexibility and approval of the gradient based iterative algorithms [28], the least squares based iterative algorithms encounter obstacles in deriving the iterative algorithms for solving the matrix equations AX + X T B = F and AXB + CXD = F [25], though it is successful at solving Ax = b, AXB = F and the coupled Sylvester matrix equations [21, 22]. How to ﬁnd a way to overcome these obstacles is the main motivation of this paper. The symmetric positive deﬁnite systems constitute one of the most important cases of special Ax = b problems [29]. To solve this kind of matrix equations, the properties of the coeﬃcient matrices have been investigated and these properties play important roles in the analysis of the convergence of the related iterative algorithms [5]. Closely related to the symmetric positive deﬁnite matrix, there is a kind of matrices and their eigenvalues have special structure and bounds which are encountered in the topics of the Jacobi iteration [6] and the preconditioning of the conjugate gradient method [5]. To our best knowledge, the properties of the eigenvalues related to this kind of matrices have not been fully investigated. This paper discusses the structure ✩ This work was supported by the National Natural Science Foundation of China (No. 61273194), the Natural Science Foundation of Jiangsu Province (China, BK2012549), and the PAPD of Jiangsu Higher Education Institutions and the 111 Project (B12018). ∗ Corresponding author Email addresses: [email protected] (Huamin Zhang), [email protected] (Feng Ding)

Preprint submitted to Journal of the Franklin Institute – Engineering and Applied Mathematics

September 10, 2013

and the bounds of the eigenvalues related to the symmetric positive deﬁnite matrices and establishes several new results to determine the range of the eigenvalues. From these results, a family of the iterative algorithms for solving linear matrix equations are obtained. Analysis shows that the least squares based iterative algorithms discussed in [21] are the special cases of this kind of the iterative algorithms. The paper is organized as follows. Section 2 gives some notations and several basic preliminary lemmas. Sections 3 establishes a new property of the eigenvalues related to the symmetric positive deﬁnitive matrices. Section 4 extends this new property to the general cases. Based on this new property, section 5 proposes a family of the iterative algorithms for solving the linear matrix equations and the relationship between this kind of iterative algorithms and the least squares based iterative algorithms [21] is discussed. Section 6 presents the least squares based iterative algorithms for solving the generalized (coupled) Sylvester matrix equations. Section 7 oﬀers an example to illustrate the eﬀectiveness of the algorithms proposed in the paper. Finally, we oﬀer some concluding remarks in Section 8. 2. Basic preliminaries Let us introduce some notation and lemmas ﬁrst. I n is the identity matrix with order n × n; I is an identity matrix of appropriate orders. For a square matrix M , λ[A] represents the set of the eigenvalues of A, λmax [A] (λmin [A]) represents the maximum (minimum) eigenvalue of A, X denotes the norm of X and is deﬁned by X2 := tr[X T X], 1m×n represents a m × n dimensional matrix whose elements are all unity For an m × n matrix X = [x1 , x2 , · · · , xn ] ∈ Rm×n , xi ∈ Rm , col[X] represents an mn-dimensional vector formed by the columns of X, i.e., col[X] := [xT1 , xT2 , · · · , xTn ]T ∈ Rmn . For two matrices A and B, A ⊗ B denotes their Kronecker product. Lemma 1. For kij ∈ R, i = 1, 2, · · · , m, j = 1, 2, · · · , n and p < n, if inequality p

0

m

j=1

i=1 n m

n m j=1 i=1

2 kij = 0, then we have the following

2 kij

j=1 i=1

2 kij

m.

Proof is easy and omitted here. The following lemma is about the block matrix determinants. Lemma 2. If A B C D

matrix A is invertible, then for any block matrix, we have = |A||D − CA−1 B|,

or if matrix D is invertible, then we have A B −1 C D = |D||A − BD C|. Lemma 3. If A ∈ Rm×n is a full column-rank matrix, then A(AT A)−1 AT is idempotent and the eigenvalues of A(AT A)−1 AT are 1 or 0, that is, there exists an orthogonal matrix Q such that QT [A(AT A)−1 AT ]Q = diag[1, · · · , 1, 0, · · · , 0] =: Λ. Furthermore, we have rank[Λ] = n. Proof. If σ ∈ λ[A(AT A)−1 AT ], then there exists a nonzero vector x ∈ Rm satisfying A(AT A)−1 AT x = σx.

2

Thus, we have [A(AT A)−1 AT x]T [A(AT A)−1 AT x] = (σx)T (σx), xT [A(AT A)−1 AT ][A(AT A)−1 AT ]x = σ 2 x2 , xT [A(AT A)−1 AT A(AT A)−1 AT ]x = σ 2 x2 , xT [A(AT A)−1 AT ]x = σ 2 x2 , xT σx = σ 2 x2 , σx2 = σ 2 x2 . Since x2 = 0, A(AT A)−1 AT has eigenvalues σ = 0 or σ = 1. Because of the symmetry of A(AT A)−1 AT , there exists a real orthogonal matrix Q := [q 1 , q 2 , · · · , q m ] ∈ Rm×m such that QT [A(AT A)−1 AT ]Q = diag[1, · · · , 1, 0, · · · , 0] = Λ. On the other hand, since (AT A)−1 AT is the left pseudo-inverse of A, we have rank[Λ] = rank[QT [A(AT A)−1 AT ]Q] = rank[A(AT A)−1 AT ] = rank[A] = n. This proves Lemma 3. If A, B ∈ Rn×n are two symmetric positive deﬁnite matrices, and λ and u are their generalized eigenvalue and the generalized eigenvector, that is, Au = λBu. In this case, we have the following lemmas [30]. Lemma 4. If A and B are symmetric positive definite matrices, then the eigenvalues of A−1 B are positive. Lemma 5. (Generalized Real Schur Decomposition) If the eigenvalues of A, B ∈ Rn×n are real, then there exist orthogonal matrices Q and Z such that QT AZ is upper triangular and QT BZ is upper triangular. 3. A property of the symmetric positive deﬁnite matrices In this section, we establish a new property about the symmetric positive deﬁnite matrix. If M ∈ Rt×t is a symmetric positive deﬁnite matrix, then there exists an invertible matrix A such that M = AAT . Suppose that A can be expressed as a block matrix A1 A11 A12 A := := , A2 A21 A22 where A11 ∈ Rm×m , A12 ∈ Rm×n , A21 ∈ Rn×m , and A22 ∈ Rn×n ; in this case, M can be expressed as A1 A1 AT1 A1 AT2 T T T AA = [A1 , A2 ] = . A2 A2 AT1 A2 AT2 Set

N :=

A1 AT1 0

0 A2 AT2

∈ Rt×t .

With these notations, we have the following conclusion. Theorem 1. If δ ∈ λ[N −1 M ], then 0 < δ < 2. Proof. Let p(x) := |xI t − N −1 M | be the characteristic polynomial of matrix N −1 M , we have Im (A1 AT1 )−1 (A1 AT2 ) p(x) = xI t − . (A2 AT2 )−1 (A2 AT1 ) In We show that 2 ∈ λ[N −1 M ]. If 2 ∈ λ[N −1 M ], then Im −(A1 AT1 )−1 (A1 AT2 ) p(2) = = 0. −(A2 AT2 )−1 (A2 AT1 ) In

3

(1)

On the other hand, Im −(A1 AT1 )−1 (A1 AT2 ) p(2) = In −(A2 AT2 )−1 (A2 AT1 ) (A1 AT1 )−1 0 −A1 AT2 A1 AT1 = 0 (A2 AT2 )−1 −A2 AT1 A2 AT2 (A1 AT1 )−1 A1 0 = [AT1 , −AT2 ] = 0. T −1 −A2 0 (A2 A2 ) This is a contradiction. Hence 2 ∈ λ[N −1 M ]. From Lemma 2, suppose that m n, Equation (1) can (x − 1)I m −(A1 AT1 )−1 (A1 AT2 ) p(x) = T −1 T (x − 1)I n −(A2 A2 ) (A2 A1 )

be rewritten as

= |(x − 1)I m ||(x − 1)I n − (A2 AT2 )−1 (A2 AT1 )(x − 1)−1 I m (A1 AT1 )−1 (A1 AT2 )| = |(x − 1)I m ||(x − 1)I n − (x − 1)−1 I n (A2 AT2 )−1 (A2 AT1 )(A1 AT1 )−1 (A1 AT2 )| = (x − 1)m−n |(x − 1)2 I n − (A2 AT2 )−1 A2 [AT1 (A1 AT1 )−1 A1 ]AT2 |,

(2)

Since AT1 (A1 AT1 )−1 A1 is a real symmetric matrix, according to Lemma 3, there exists an orthogonal matrix Q := [q 1 , q 2 , · · · , q t ] such that [AT1 (A1 AT1 )−1 A1 ][q 1 , q 2 , · · · , qt ] = [q 1 , q 2 , · · · , q t ]diag[1, · · · , 1, 0, · · · , 0] = [q 1 , q 2 , · · · , q m , 0, 0, · · · , 0],

(3)

where m = rank[AT1 (A1 AT1 )−1 A1 ]. Let P := (A2 AT2 )−1 A2 [AT1 (A1 AT1 )−1 A1 ]AT2 ∈ Rn×n and suppose that ρ ∈ λ[P ]. There exists a nonzero vector y ∈ Rn such that P y = ρy, this equation can be rewritten as (A2 AT2 )−1 A2 [AT1 (A1 AT1 )−1 A1 ]AT2 y = ρy, A2 [AT1 (A1 AT1 )−1 A1 ]AT2 y = ρ(A2 AT2 )y.

(4)

It is clear that AT2 y can be expressed as a linear combination of q 1 , q 2 , · · · , q t , i.e., AT2 y = k1 q 1 + k2 q 2 + · · · + kt q t ,

(5)

where ki ∈ R, i = 1, 2, · · · , t and k12 + k22 + · · · + kt2 = 0. Using (3) and (5), we have (AT2 y)T (AT2 y) = (k1 q 1 + k2 q 2 + · · · + kt q t )T (k1 q 1 + k2 q 2 + · · · + kt q t ) = k12 + k22 + · · · + kt2 , −1

[A1 (A1 A1 ) T

T

T

T

(6)

−1

A1 ]A2 y = [A1 (A1 A1 ) A1 ](k1 q 1 + k2 q 2 + · · · + kt q t ) = k1 q 1 + k2 q 2 + · · · + km q m . T

(7)

Together with (4), (6) and (7), it gives y T A2 [AT1 (A1 AT1 )−1 A1 ]AT2 y = ρy T A2 AT2 y,

y T A2 (k1 q 1 + k2 q 2 + · · · + km q m ) = ρ(k12 + k22 + · · · + kt2 ), (k1 q 1 + k2 q 2 + · · · + kt q t )T (k1 q 1 + k2 q 2 + · · · + km q m ) = ρ(k12 + k22 + · · · + kt2 ), 2 (k12 + k22 + · · · + km ) = ρ(k12 + k22 + · · · + kt2 ).

Hence, we have 0ρ=

2 k12 + k22 + · · · + km 1. 2 2 k1 + k2 + · · · + kt2

(8)

Since 0 ρ 1, we suppose that λ[P ] = {ρ21 , · · · , ρ2n } with 0 ρn · · · ρ1 1. According to the Schur decomposition theorem, there exists an orthogonal matrix U such that U T P U = diag[ρ21 , · · · , ρ2n ] + R,

(9)

4

where R is strictly upper triangular. Substituting (9) into (2) and simplifying it give p(x) = (x − 1)m−n |(x − 1)2 I n − P |, = (x − 1)m−n |U T ||(x − 1)2 I n − P ||U | = (x − 1)m−n |U T (x − 1)2 I n U − U T P U | = (x − 1)m−n (x − 1)2 I n − diag[ρ21 , · · · , ρ2n ] − R = (x − 1)m−n [(x − 1)2 − ρ21 ] · · · [(x − 1)2 − ρ2n ] = (x − 1)m−n (x − 1 − ρ1 )(x − 1 + ρ1 ) · · · (x − 1 − ρn )(x − 1 + ρn ). From this equation, we can see that the eigenvalues of N −1 M are 1 + ρ1 , 1 − ρ1 , · · · , 1 + ρn , 1 − ρn , 1, · · · , 1.

(10)

For 2 ∈ λ[N −1 M ], an improvement of inequality (8) is 0 ρ < 1. From this inequality, we have 0 < δ < 2. This proves theorem 1. From (10), we have λmax [N −1 M ] + λmin [N −1 M ] = 2.

(11)

This result will be used in Section 5. From the above discussion, we can ﬁnd that the supposition m n is not essential, in fact, if m < n then there exists the permutation matrix S such that (x − 1)I m −(A1 AT1 )−1 (A1 AT2 ) S T p(x) = S −(A2 AT2 )−1 (A2 AT1 ) (x − 1)I n (x − 1)I n −(A2 AT2 )−1 (A2 AT1 ) = −(A1 AT1 )−1 (A1 AT2 ) (x − 1)I m = |(x − 1)I n ||(x − 1)I m − (A1 AT1 )−1 (A1 AT2 )(x − 1)−1 I n (A2 AT2 )−1 (A2 AT1 )| = |(x − 1)I n ||(x − 1)I m − (x − 1)−1 I m (A1 AT1 )−1 A1 [AT2 (A2 AT2 )−1 A2 ]AT1 | = (x − 1)n−m (x − 1)2 I m − (A1 AT1 )−1 A1 [AT2 (A2 AT2 )−1 A2 ]AT1 . This transformation does not change the subsequent proof. If M ∈ Rt×t is a symmetric positive deﬁnite matrix, M (m) ∈ Rm×m is a principal submatrix of M and M (m) ∈ R(t−m)×(t−m) is the complementary submatrix of M (m), then there exists a permutation matrix S such that M (m) ∗ , SM S T = ∗ M (m) where ∗ denote the proper block matrices. We have the following theorem. Theorem 2. If M ∈ Rt×t is a symmetric positive definite matrix, let M (m) 0 N := , 0 M (m) where M (m) ∈ Rm×m , 1 m t, is a leading principle submatrix of M and M (m) ∈ R(t−m)×(t−m) is the complementary submatrix of M (m). In this case, if δ ∈ λ[N −1 M ], then we have 0 < δ < 2. As a special case, if m = t, then N = M , N −1 M = I t and δ = 1. Remark 1: Compared with the classical Gerschgorin circle theorem in [6], the conclusions of Theorems 1 and 2 are very useful and interesting. For example, consider the eigenvalues λ1 and λ2 of matrix A=

1 −0.5 −0.25 1

=

2 0 0 4

−1

2 −1 −1 4

.

According to the Gerschgorin circle theorem, we can only determine the eigenvalue λ2 within the bounds of 0.5 < λ2 < 1.5 even if we know that λ1 = 0.6464. But from Theorem 1, we can easily get λ2 = 1.3536 because of the fact of λ1 + λ2 = 2 from Theorem 1.

5

4. General cases Next, we study the general cases of Theorem 1. If M ∈ Rt×t is a symmetric positive deﬁnite matrix, then there exists an invertible matrix A such that M = AAT . Suppose that A can be expressed as a block matrix, ⎤ ⎡ ⎤ ⎡ A11 A12 A13 A1 A := ⎣ A2 ⎦ := ⎣ A21 A22 A23 ⎦ , A3 A31 A32 A33 where A11 ∈ Rm×m , A22 ∈ Rn×n and A33 ∈ Rp×p are proper block matrices. M can be expressed as ⎤ ⎡ ⎡ A1 A1 AT1 A1 AT2 T T T T ⎦ ⎣ ⎣ AA = A2 [A1 , A2 , A3 ] = A2 AT1 A2 AT2 A3 A3 AT1 A3 AT2 Set

⎡

A1 AT1 ⎣ 0 N := 0

0 A2 AT2 0

block square matrices, and Aij , i, j = 1, 2, 3, i = j are ⎤ A1 AT3 A2 AT3 ⎦ . A3 AT3

⎤ 0 ⎦ ∈ Rt×t . 0 T A3 A3

With these notations, we have the following theorem. Theorem 3. If δ ∈ λ[N −1 M ], then 0 < δ < 3. Proof. Let p(x) := |xI t − N −1 M |, we have ⎤ ⎡ (A1 AT1 )−1 A1 AT2 (A1 AT1 )−1 A1 AT3 Im In (A2 AT2 )−1 A2 AT3 ⎦ p(x) = xI t − ⎣ (A2 AT2 )−1 A2 AT1 T −1 T T −1 T (A3 A3 ) A3 A1 (A3 A3 ) A3 A2 Ip ⎡ ⎤ (x − 1)I m −(A1 AT1 )−1 A1 AT2 −(A1 AT1 )−1 A1 AT3 T T −1 (x − 1)I n −(A2 AT2 )−1 A2 AT3 ⎦ . = ⎣ −(A2 A2 ) A2 A1 −(A3 AT3 )−1 A3 AT1 −(A3 AT3 )−1 A3 AT2 (x − 1)I p Deﬁne P :=

(z − 1)I m −(A2 AT2 )−1 A2 AT1

−(A1 AT1 )−1 A1 AT2 (z − 1)I n

(12)

∈ R(m+n)×(m+n) ,

Q := AT3 (A3 AT3 )−1 A3 ∈ Rt×t , (A1 AT1 )−1 A1 R := Q[AT1 , AT2 ] ∈ R(m+n)×(m+n) . (A2 AT2 )−1 A2 First, we discuss the eigenvalues of R. If σ ∈ λ[R], then there exists a nonzero vector x = satisfying

x1 x2

∈ R(m+n)

Rx = σx, where x1 ∈ Rm and x2 ∈ Rn . From this equation, we have T A1 AT1 A1 A1 0 x=σ Q x, A2 A2 0 A2 AT2 T A1 A1 A1 AT1 0 T T x Q x, x = σx 0 A2 AT2 A2 A2 T T T x1 0 A1 A1 x1 x1 A1 AT1 x1 Q =σ , 0 A2 AT2 x2 A2 A2 x2 x2 x2 (AT1 x1 + AT2 x2 )T Q(AT1 x1 + AT2 x2 ) = σ[(AT1 x1 )T (AT1 x1 ) + (AT2 x2 )T (AT2 x2 )]. Since Q is idempotent, there exists a real orthogonal S := [s1 , s2 , · · · , st ] such that S T QS = diag[1, · · · , 1, 0, · · · , 0].

6

(13)

It is clear that AT1 x1 and AT2 x2 can be expressed as a linear combination of s1 , s2 , · · · , st , i.e., AT1 x1 = k1 s1 + k2 s2 + · · · + kt st , AT2 x2 = l1 s1 + l2 s2 + · · · + lt st . Adding these two equations together, we have AT1 x1 + AT2 x2 = (k1 + l1 )s1 + (k2 + l2 )s2 + · · · + (kt + lt )st , (A1 x1 ) (A1 x1 ) + (A2 x2 ) T

T

T

T

T

(A2 x2 ) = (k12 T

+

k22

+ ···+

kt2 )

+

(l12

+

l22

+ ··· +

(14)

lt2 ).

(15)

By Lemma 3, the eigenvalues of Q are 1,0 and rank[Q] = p. Substituting (14), (15) into (13) gives (k1 + l1 )2 + (k2 + l2 )2 + · · · + (kp + lp )2 = σ(k12 + l12 + k22 + l22 + · · · + kt2 + lt2 ). From Lemma 1, we have 0σ=

(k1 + l1 )2 + (k2 + l2 )2 + · · · + (kp + lp )2 2. k12 + l12 + k22 + l22 + · · · + kt2 + lt2

That is, if σ ∈ λ[R], then 0 σ 2. Second, we study the characteristic polynomial p(x) of matrix N −1 M . To simplify p(x), we suppose that p m + n. From Lemma 2, we have −(A1 AT1 )−1 A1 AT3 T −1 T T −1 T −1 p(z) = |(z − 1)I p | P − I [−(A A ) A A , −(A A ) A A ] (z − 1) T −1 T p 3 3 3 1 3 3 3 2 −(A2 A2 ) A2 A3 = |(z − 1)I p | P − (z − 1)−1 I m+n R = (z − 1)p−(m+n) |(z − 1)I m+n P − R| .

(16)

Since N and M are symmetric positive deﬁnite matrices, from Lemma 4, we can see that the eigenvalues of N −1 M are positive. From Theorem 1 and Lemma 5, Equation (16) can be rewritten as p(x) = (x − 1)p−(m+n) |(x − 1)I m+n P − R]| = (x − 1)p−(m+n) |U V T | ⎡ (x − 1)(x − ρ1 ) ∗ ··· ⎢ 0 (x − 1)(x − ρ ) ··· 2 ⎢ × ⎢ .. .. .. ⎣ . . . 0 0 ··· p−(m+n)

= ±(x − 1)

∗ ∗ .. .

(x − 1)(x − ρm+n )

⎤

⎡

⎥ ⎢ ⎥ ⎢ ⎥−⎢ ⎦ ⎣

σ1 0 .. .

∗ σ2 .. .

0

0

··· ··· .. .

···

∗ ∗ .. . σm+n

⎤ ⎥ ⎥ ⎥ ⎦

[(x − 1)(x − ρ1 ) − σ1 ][(x − 1)(x − ρ2 ) − σ2 ] · · · [(x − 1)(x − ρm+n ) − σm+n ],

where U and V are orthogonal matrices and ρi ∈ λ[P ] and σi ∈ λ[R], i = 1, 2, · · · , m + n. By studying the equation (z − 1)(z − ρi ) − σi = 0, i = 1, 2, · · · , m + n, we have 0 < δ < 3. On the other hand, if p < m + n, then T 0 (A1 AT1 )−1 A1 A1 rank[R] = rank Q 0 (A2 AT2 )−1 A2 A2 T A1 A = rank Q 1 A2 A2 T A1 A1 = rank AT3 (A3 AT3 )−1 A3 A2 A2 T A1 A1 = rank AT3 A3 A2 A2 T A1 = rank A3 A2 ≤ p.

(17)

Simplifying (12) gives

−(A1 AT1 )−1 A1 AT3 T −1 T T −1 T −1 p(x) = |(x − 1)I p | P − I [−(A A ) A A , −(A A ) A A ] (x − 1) p 3 3 3 1 3 3 3 2 T −1 T −(A2 A2 ) A2 A3 = |(x − 1)I p | P − (x − 1)−1 I m+n R . 7

From Theorem 1, Lemmas 4 and 5, we have −(A1 AT1 )−1 A1 AT3 T −1 T T −1 T −1 p(x) = |(x − 1)I p | P − (x − 1) I p [−(A3 A3 ) A3 A1 , −(A3 A3 ) A3 A2 ] T −1 T −(A2 A2 ) A2 A3 = |(x − 1)I p | P − (x − 1)−1 I m+n R = |(x − 1)I p ||U V T | × diag[(x − ρ1 ), (x − ρ2 ), · · · , (x − ρm+n )] + W − (x − 1)−1 I m+n diag[σ1 , · · · , σp , 0, · · · , 0] − Y = ± diag[(x − 1), · · · , (x − 1), 1, · · · , 1] p m+n × diag[(x − ρ1 ), (x − ρ2 ), · · · , (x − ρm+n )] − (x − 1)−1 I m+n diag[σ1 , · · · , σp , 0, · · · , 0] = ± diag[(x − 1)(x − ρ1 ) · · · (x − 1)(x − ρp )(x − ρp+1 ), · · · , (x − ρm+n )] − diag[σ1 , · · · , σp , 0, · · · , 0] p = ±[(x − 1)(x − ρ1 ) − σ1 ] · · · [(x − 1)(x − ρp ) − σp ](x − ρp+1 ) · · · (x − ρm+n ), where U and V are orthogonal matrices, W and Y are two proper strictly upper triangular matrices. By analyzing equation (z − 1)(z − ρi ) − σi = 0, i = 1, 2, · · · , p and x − ρj = 0, j = p + 1, · · · , m + n, we have 0 < δ < 3. This proves Theorem 3. From the above discussion, we have the next theorem expressing a more general situation. Theorem 4. If M ∈ Rt×t is a symmetric positive definite matrix, and M (1), M (2), · · · , M (n) are the principle submatrices, where M (i) ∈ Rti ×ti , t1 + t2 + · · · + tn = t, then there exists a permutation matrix S such that ⎡ ⎤ M (1) ∗ ··· ∗ ⎢ ⎥ ∗ M (2) · · · ∗ ⎢ ⎥ S TM S = ⎢ ⎥, .. .. .. .. ⎣ ⎦ . . . . ∗ ∗ · · · M (n) where ∗ denote the proper submatrices. In ⎡ M (1) 0 ··· 0 ⎢ 0 M (2) · · · 0 ⎢ N := ⎢ .. .. .. .. ⎣ . . . . 0

0

···

this case, set ⎤ ⎥ ⎥ ⎥, ⎦

M (n)

if δ ∈ λ[N −1 M ], then 0 < δ < n. Proof is easy and omitted here. Remark 2: The structure and range of this kind of matrix is often encountered in the numerical linear algebra [5, 6]. The conjugate gradient method is eﬀective for solving large-scale linear systems. In order to improve the rate of convergence, the concept of the precondition matrix was introduced [5]. Consider the linear system M x = b, where M is a symmetric positive deﬁnite matrix. If we pre-multiply this system by N −1 , then the coeﬃcient matrix of this new system is N −1 M . From Theorem 4, the proper choice of N −1 will contribute to the desired range and distribution of the eigenvalues of N −1 M . This precondition transformation will improve the eﬃciency of the conjugate gradient method obviously [5]. So N −1 is one of the important precondition matrices and Theorems 1 to 4 play important roles in the numerical linear algebra. 5. The iterative algorithms for linear matrix equations From Theorems 2 to 4, we can obtain a family of the iterative algorithms for linear matrix equations. Consider the matrix equation AX = F , where A ∈ Rm×n , F ∈ Rm×p are constant matrices and X ∈ Rn×p is unknown.

8

5.1. The matrix equation AX = F T n×n Theorem 5. For the equation and AX = F , if A is a full column-rank matrix, let M := A A ∈ R M (m) 0 N := , where M (m) ∈ Rm×m , 1 ≤ m < n is a leading principle submatrix of M and 0 M (n − m) M (n − m) ∈ R(n−m)×(n−m) is the complementary submatrix of M (m), then the iterative algorithm X(k) = X(k − 1) + μN −1 AT [F − AX(k − 1)], 2 0 < μ< , λmax [N −1 M ]

(18)

yields limk→∞ X(k) = X. In this case, the best convergence factor is 1. ˜ Proof. Deﬁning the estimation error matrix X(k) := X(k) − X, substituting F = AX into (18) and using ˜ X(k) = X(k) − X, we obtain ˜ ˜ − 1) − μN −1 M X(k ˜ − 1) X(k) = X(k ˜ − 1). = (I n − μN −1 M )X(k From Theorem 2, if δi ∈ λ[N −1 M ], then 0 < δi < 2. The eigenvalues of I n − μN −1 M are 1 − μδi . If ˜ = 0. In this case, we have 0 < μ < δ2i , and their intersection is −1 < 1 − μδi < 1, then limk→∞ X(k) 2 . The best factor μ0 satisﬁes 0<μ< −1 λmax [N M ] |1 − μ0 λmax [N −1 M ]| = |1 − μ0 λmin [N −1 M ]|, From (10) and (11), solving this equation gives the best step-size μ0 =

2 = 1. λmax [N −1 M ] + λmin [N −1 M ]

This proves Theorem 5. Remark 3: Some interpretations will expose the relationship between the iterative algorithm in (18) and the Jacobi iteration. Consider the matrix equation Ax = b, where A and b are known, and A = [aij ] ∈ Rm×m is a symmetric positive deﬁnite matrix, x is the vector to be determined. Under this assumption, a simple version of algorithm (18) is x(k) = x(k − 1) + μN −1 [b − Ax(k − 1)], 2 . 0 < μ< λmax [N −1 A]

(19)

If we take μ = 1 and N = diag[a11 , a22 , · · · , ann ], then we quickly ﬁnd that algorithm (19) is the classical Jacobi iteration. Remark 4: If N is a matrix formed by the principle submatrices of A, then the iterative algorithm in (19) can be viewed as an extension of the Jacobi iteration. A simple operation shows that the performance of this iteration is superior to the Jacobi iteration. Consider the Poisson’s equation in one dimension, −

d2 v(x) = f (x), 0 < x < 1. dx2

If we split the interval (0, 1) into n small sections and after some manipulations, we transform this problem to solve a matrix equation Ax = b [7]. Here, b is known and A ∈ Rn×n is a symmetric tridiagonal matrix ⎡ ⎤ 2 −1 0 ⎢ ⎥ ⎢ −1 . . . . . . ⎥ ⎥. A=⎢ ⎢ ⎥ .. .. ⎣ . . −1 ⎦ 0 −1 2 For simplicity, let n = 5 and N 1 = diag[M 1 , M 2 ], where ⎡ ⎤ 2 −1 0 2 −1 ⎣ ⎦ M 1 = −1 2 −1 , M 2 = . −1 2 0 −1 2 Let ρ[A] := max{|λ| : λ ∈ λ[A]} denotes the spectral radius of A. Taking μ = 1, we have ρ[I 5 −N −1 1 A] = 0.7071. But according to the Jacobi iteration, we have N 2 = diag[2, 2, · · · , 2] and obtain ρ[I 5 − N −1 2 A] = 0.8660. Since −1 ρ[I 5 − N −1 1 A] < ρ[I 5 − N 2 A], the iterative algorithm in (19) is an improvement and extension of the Jacobi iteration. 9

5.2. The coupled Sylvester matrix equations In this subsection, we extend the algorithm in Theorem 5 to the coupled Sylvester matrix equations, AX + Y B = C, DX + Y E = F ,

(20)

where A, D ∈ Rm×m , B, E ∈ Rn×n and C, F ∈ Rm×n are given constant matrices, X, Y ∈ Rm×n are the unknown matrices to be solved. From Theorem 5 and using the Kronecker product, if equation (20) has the unique solution, then the corresponding iterative algorithm with best convergence factor μ = 1 can be expressed as T A C − Y (k − 1)B , (21) X(k) = (AT A + DT D)−1 D F − Y (k − 1)E −1

Y (k) = [C − AX(k − 1), F − DX(k − 1)][B, E]T (BB T + EE T )

.

(22)

The least squares iterative algorithm for solving matrix equations equation (20), which was presented in [21], can be expressed as T C − Y (k − 1)B − AX(k − 1) T T −1 A X(k) = X(k − 1) + μ(A A + D D) , (23) D F − Y (k − 1)E − DX(k − 1) Y (k) = Y (k − 1) + μ[C − AX(k − 1) − Y (k − 1)B, F − DX(k − 1) − Y (k − 1)E] −1

×[B, E]T (BB T + EE T ) 2 0<μ< . λmax [N −1 M ]

,

(24) (25)

Here, N = diag[I n ⊗ (AT A + D T D), (BB T + EE T ) ⊗ I m ] and I n ⊗ (AT A + D T D) B T ⊗ AT + E T ⊗ DT M= . B⊗A+E⊗D (BB T + EE T ) ⊗ I m It is not diﬃcult to ﬁnd that if μ = 1 then algorithm (23)–(25) is equivalent to the algorithm in (21)–(22). From the above discussion, we can ﬁnd that if matrix equations or coupled matrix equations [21] have an unique solution then it always can be written together as an equation by using the Kronecker product, in this case, the algorithm in Theorem 5 is always eﬀective. 6. The general coupled matrix equations Consider the following coupled matrix equations, A1 XB 1 + C 1 Y D 1 = F 1 , A2 XB 2 + C 2 Y D 2 = F 2 .

(26)

By the formula col[AXB] = (B T ⊗ A)col[X], system (26) can be expressed as T B 1 ⊗ A1 DT1 ⊗ C 1 col[X, Y ] = col[F 1 , F 2 ]. B T2 ⊗ A2 DT2 ⊗ C 2 Set

A :=

B T1 ⊗ A1 B T2 ⊗ A2

and thus we have

DT1 ⊗ C 1 DT2 ⊗ C 2

B 1 ⊗ AT1 D 1 ⊗ C T1

(27)

,

B 2 ⊗ AT2 D2 ⊗ C T2

B T1 ⊗ A1 B T2 ⊗ A2

DT1 ⊗ C 1 DT2 ⊗ C 2

, M := A A = (B 1 B T1 ) ⊗ (AT1 A1 ) + (B 2 B T2 ) ⊗ (AT2 A2 ) (B 1 DT1 ) ⊗ (AT1 C 1 ) + (B 2 D T2 ) ⊗ (AT2 C 2 ) = . (D 1 B T1 ) ⊗ (C T1 A1 ) + (D 2 B T2 ) ⊗ (C T2 A2 ) (D 1 D T1 ) ⊗ (C T1 C 1 ) + (D 2 DT2 ) ⊗ (C T2 C 2 ) T

Let

N :=

(B 1 B T1 ) ⊗ (AT1 A1 ) + (B 2 B T2 ) ⊗ (AT2 A2 ) 0 0 (D 1 D T1 ) ⊗ (C T1 C 1 ) + (D 2 DT2 ) ⊗ (C T2 C 2 ) 10

.

According to Theorem 5, we have the following iterative algorithm for system (27), col[X(k), Y (k)] = col[X(k − 1), Y (k − 1)] +μN −1 AT {col[F 1 , F 2 ] − Acol[X(k − 1), Y (k − 1)]}, 2 0<μ< , λmax [N −1 M ]

(28)

It is clear that from the iterative solutions of equation (27), we can easily get the solution for equation (26) . Consider the generalized Sylvester matrix equation AXB + CXD = F , using the formula col[AXB] = (B T ⊗ A)col[X], this equation can be rewritten as (B T ⊗ A + DT ⊗ C)col[X] = col[F ].

(29)

If B, D ∈ R2n×2n , then these two matrices can be expressed as D11 D12 B 11 B 12 , D= , B= B 21 B 22 D21 D22 where B ij , D ij ∈ Rn×n . Let M := (B ⊗ AT + D ⊗ C T )(B T ⊗ A + D T ⊗ C) = (BB T ) ⊗ (AT A) + (DB T ) ⊗ (C T A) + (BD T ) ⊗ (AT C) + (DD T ) ⊗ (C T C), and

0 B 11 B T11 + B 12 B T12 ⊗ (AT A) 0 B 21 B T21 + B 22 B T22 D11 B T11 + D 12 B T12 0 + ⊗ (C T A) 0 D 21 B T21 + D22 B T22 0 B 11 DT11 + B 12 D T12 + ⊗ (AT C) 0 B 21 DT21 + B 22 D T22 D11 D T11 + D12 DT12 0 + ⊗ (C T C). 0 D21 D T21 + D22 DT22

N :=

According to Theorem 5, we can get the following iterative algorithm for solving equation (29): col[X(k)] = col[X(k − 1)] + μN −1 (B ⊗ AT + D ⊗ C T ) × {col[F ] − (B T ⊗ A + D T ⊗ C)col[X(k − 1)]} , 2 0 < μ< . λmax [N −1 M ]

(30) (31)

Thus, we can easily get the solution for equation AXB + CXD = F from the iterative solution of equation (29). 7. Example In this section, we oﬀer an example to illustrate the performance of the iterative algorithm in this paper. Suppose that the coupled matrix equations are AX + Y B = C, DX + Y E = F with 2 7 −10 0 10 −47 A= , B= , C= , −11 28 −6 3 4 −55 −2 7 −6 24 −178 13 D= , E= , F = . 1 −3 20 17 80 230 The exact solution is x11 x12 2 −1 X= = , x21 x22 4 −3

Y =

y11 y21

y12 y22

=

7 5

−8 6

.

From (25), we have λmax [N −1 M ] = 1.4779 and 0 < μ < 2/(λmax [N −1 M ])−1 = 1.3533. Taking X(0) = 10−6 12×2 , we apply the iterative algorithm in (23)–(24) to compute X(k) and Y (k), the iterative solutions of X and Y are shown in Tables 1–3 with the relative errors X(k) − X2 + Y (k) − Y 2 × 100%. δ := X2 + Y 2 11

Table 1: The least squares based iterative solutions with μ = 0.7

k 1 2 3 4 5 10 15 20 Solution

x11 -1.40857 -1.06180 -0.34654 0.33730 0.86145 1.85556 1.98325 1.99813 2.00000

k 1 2 3 4 5 10 15 20 Solution

x11 -2.01224 -0.44217 0.66696 1.38237 1.62577 1.99189 1.99964 1.99999 2.00000

k 1 2 3 4 5 10 15 20 Solution

x11 -2.41469 0.32915 1.24477 1.86893 1.67158 2.19336 1.91828 2.02817 2.00000

x12 -4.77985 -3.36731 -2.39167 -1.81008 -1.47255 -1.03257 -1.00210 -1.00011 -1.00000

x21 -1.51612 0.32108 1.57831 2.43326 2.99452 3.89474 3.98913 3.99888 4.00000

x22 -3.60733 -3.54018 -3.34862 -3.20687 -3.11801 -3.00445 -2.99977 -2.99992 -3.00000

y11 3.78579 4.90082 5.65453 6.13878 6.45027 6.94287 6.99411 6.99939 7.00000

y12 -5.17092 -6.28654 -6.94572 -7.34743 -7.59612 -7.96380 -7.99675 -7.99971 -8.00000

y21 3.17009 3.72510 4.17581 4.48281 4.68035 4.97317 4.99780 4.99982 5.00000

y22 3.57154 4.65721 5.23865 5.55520 5.73440 5.97669 5.99781 5.99979 6.00000

δ (%) 64.27687 44.04547 29.38079 19.23939 12.49031 1.38019 0.14983 0.01615

y22 5.10220 5.13121 5.86434 5.82253 5.97123 5.99832 5.99998 6.00000 6.00000

δ (%) 68.83497 27.59522 15.58887 6.37102 3.74090 0.07545 0.00295 0.00005

y22 6.12263 4.93988 6.43494 5.54415 6.29603 5.92830 6.01683 5.99594 6.00000

δ (%) 78.63425 35.61169 26.00183 16.71051 12.65214 3.04022 0.87133 0.25952

Table 2: The least squares based iterative solutions with μ = 1.0

x12 -6.82835 -1.01919 -1.89985 -1.00875 -1.15632 -0.99966 -1.00001 -1.00000 -1.00000

x21 -2.16589 2.51174 2.56554 3.65023 3.67115 3.99581 3.99979 4.00000 4.00000

x22 -5.15333 -2.80772 -3.31323 -2.95973 -3.04461 -2.99927 -2.99998 -3.00000 -3.00000

y11 5.40827 5.36601 6.60628 6.62610 6.90944 6.99553 6.99994 7.00000 7.00000

y12 -7.38703 -6.49793 -7.82505 -7.70023 -7.96214 -7.99748 -7.99998 -8.00000 -8.00000

y21 4.52870 3.72049 4.89547 4.75226 4.97994 4.99826 4.99999 5.00000 5.00000

Table 3: The least squares based iterative solutions with μ = 1.2

x12 -8.19403 1.80997 -3.38563 0.44457 -2.04686 -0.81906 -1.03108 -0.99516 -1.00000

x21 -2.59907 4.65653 1.95123 4.97374 3.02902 4.24866 3.93058 4.01931 4.00000

x22 -6.18400 -1.56952 -4.04147 -2.32955 -3.45480 -2.93128 -3.00835 -2.99980 -3.00000

y11 6.48992 5.13109 7.55820 6.25419 7.45303 6.86379 7.03773 6.98953 7.00000

y12 -8.86443 -5.81124 -8.97120 -7.11188 -8.57617 -7.87347 -8.02765 -7.99380 -8.00000

y21 5.43444 3.18373 5.89700 4.19870 5.53601 4.88968 5.02202 4.99559 5.00000

From Figure 7 and Tables 1–3, we can see that the relative errors δ tends to zero as the iteration time k increases and that the best convergence factor μ equals 1. In the simulation, we keep changing the convergence factor μ from 0.1 to 1.35 and ﬁnd that the iterative solutions always converge to the exact ones. If μ ∈ (0, 1.3533) then the iterative solutions will diverge. This shows that the range of μ is appropriate. 8. Conclusions This paper discusses the structure and the bounds of the eigenvalues which related to the symmetric positive deﬁnite matrices. Several theorems have been established. These results can be used in the estimation of the eigenvalues of the related matrices. Based on these results, we present a family of iterative algorithms for solving linear matrix equations and coupled matrix equations. These iterative algorithms include the Jacobi iteration and the least squares iterative algorithms as their special cases. The analysis indicates that the iterative solutions given by the proposed algorithms can converge fast to their exact solutions under proper conditions. An example demonstrates that the algorithm is eﬀective and the range of the convergence factor is appropriate. The iterative algorithms are very related recursive or iterative parameter estimation methods in system identiﬁcation, including the maximum likelihood identiﬁcation methods [31, 32, 33], the hierarchical identiﬁcation methods [34, 35], the iterative identiﬁcation methods [36, 37, 38], and other identiﬁcation methods [39, 40, 41, 42, 43]. References [1] F. Ding, System Identiﬁcation – New Theory and Methods, Beijing: Science Press, 2013.

12

1

δ

0.8

0.6

0.4

μ = 0.70

0.2

0

μ = 1.20 μ = 1.00 0

2

4

6

8

10

12

14

16

18

20

k

Figure 1: The errors δ versus k

[2] F. Ding, Decomposition based fast least squares algorithm for output error systems, Signal Processing 93 (5) (2013) 1235-1242. [3] F. Ding, Coupled-least-squares identiﬁcation for multivariable systems, IET Control Theory and Applications 7 (1) (2013) 68-79. [4] J.Z. Liu, J. Zhang, Y. Liu, New solution bounds for the continuous algebraic Riccati equation, Journal of the Franklin Institute 348 (8)(2011) 2128-2141. [5] G.H. Golub, C.F. Van Loan, Matrix Computations, 3rd ed. Baltimore, MD: Johns Hopkins University Press, 1996. [6] R.S. Varga, Matrix Iterative Analysis, 2nd ed. Springer, 2000. [7] J.W. Demmel, Applied Numerical Linear Algebra, Society for Industrial Mathematics, 1997. [8] K.F. Liang, J.Z. Liu, Iterative algorithms for the minimum-norm solution and the least-squares solution of the linear matrix equations A1 XB1 + C1 X T D1 = M1 , A2 XB2 + C2 X T D2 = M2 , Applied Mathematics and Computation 218 (7) (2011) 3166-3175. [9] M. Dehghan, M. Hajarian, Two algorithms for ﬁnding the Hermitian reﬂexive and skew-Hermitian solutions of Sylvester matrix equations, Applied Mathematics Letters 24 (4) (2011) 444-449. [10] M. Dehghan, M. Hajarian, Iterative algorithms for the generalized centro-symmetric and central anti-symmetric solutions of general coupled matrix equations, Engineering Computations 29 (5) (2012) 528-560. [11] F. Ding, X.G. Liu, J. Chu, Gradient-based and least-squares-based iterative algorithms for Hammerstein systems using the hierarchical identiﬁcation principle, IET Control Theory and Applications 7 (2) (2013) 176-184. [12] F. Ding, H.H. Duan, Two-stage parameter estimation algorithms for Box-Jenkins systems, IET Signal Processing 7 (2013). doi: 10.1049/iet-spr.2012.0183 [13] F. Ding, Hierarchical multi-innovation stochastic gradient algorithm for Hammerstein nonlinear system modeling, Applied Mathematical Modelling 37 (4) (2013) 1694-1704. [14] F. Ding, Two-stage least squares based iterative estimation algorithm for CARARMA system modeling. Applied Mathematical Modelling 37 (7) (2013) 4798-4808. [15] J. Ding, L.L. Han, X.M. Chen, Time series AR modeling with missing observations based on the polynomial transformation, Mathematical and Computer Modelling 51 (5-6) (2010) 527-536. [16] Y.J. Liu, Y.S. Xiao, X.L. Zhao, Multi-innovation stochastic gradient algorithm for multiple-input single-output systems using the auxiliary model, Applied Mathematics and Computation 215 (4) (2009) 1477-1483. [17] Y.J. Liu, J. Sheng, R.F. Ding, Convergence of stochastic gradient estimation algorithm for multivariable ARX-like systems, Computers & Mathematics with Applications 59 (8) (2010) 2615-2627. [18] F. Ding, Combined state and least squares parameter estimation algorithms for dynamic systems. Applied Mathematical Modelling, 2013, 37(x): doi: 10.1016/j.apm.2013.06.007 [19] J. Ding, F. Ding, Bias compensation based parameter estimation for output error moving average systems, International Journal of Adaptive Control and Signal Processing 25 (12) (2011) 1100-1111. [20] F. Ding, T. Chen, Gradient based iterative algorithms for solving a class of matrix equations, IEEE Transactions on Automatic Control 50 (8) (2005) 1216-1221. [21] F. Ding, T. Chen, Iterative least squares solutions of coupled Sylvester matrix equations, Systems & Control Letters 54 (2) (2005) 95-107. [22] F. Ding, T. Chen, On iterative solutions of general coupled matrix equations, SIAM Journal on Control and Optimization 44 (6) (2006) 2269-2284. [23] L. Xie, J. Ding, F. Ding, Gradient based iterative solutions for general linear matrix equations, Computers & Mathematics with Applications 58 (7) (2009) 1441-1448. [24] J. Ding, Y.J. Liu, F. Ding, Iterative solutions to matrix equations of form Ai XBi = Fi , Computers & Mathematics with Applications 59 (11) (2010) 3500-3507. [25] F. Ding, X.P. Liu, J. Ding, Iterative solutions of the generalized Sylvester matrix equations by using the hierarchical identiﬁcation principle, Applied Mathematics and Computation 197 (1) (2008) 41-50.

13

[26] A.G. Wu, L.L. Lv, G.R. Duan, Iterative algorithms for solving a class of complex conjugate and transpose matrix equations, Applied Mathematics and Computation 217 (21) (2011) 8343-8353. [27] C.Q. Song, G.L. Chen, L.L. Zhao, Iterative solutions to coupled Sylvester-transpose matrix equations, Applied Mathematical Modelling 35 (10) (2011) 4675-4683. [28] M. Dehghan, M. Hajarian, Analysis of an iterative algorithm to solve the generalized coupled Sylvester matrix equations, Applied Mathematical Modelling 35 (7) (2011) 3285-3300. [29] X.Z. Zhong, T. Zhang, Y. Shi, Oscillation and nonoscillation of neutral diﬀerence equation with positive and negative coeﬃcients, International Journal of Innovative Computing, Information and Control 5(5) (2009) 1329-1342. [30] M.C. Pease, Methods of Matrix Algebra, New York, Academic Press, 1965. [31] J.H. Li, F. Ding, G.W. Yang, Maximum likelihood least squares identiﬁcation method for input nonlinear ﬁnite impulse response moving average systems, Mathematical and Computer Modelling 55 (3-4) (2012) 442-450. [32] J.H. Li, F. Ding, Maximum likelihood stochastic gradient estimation for Hammerstein systems with colored noise based on the key term separation technique, Computers & Mathematics with Applications 62 (11) (2011) 4170-4177. [33] W. Wang, F. Ding, J.Y. Dai, Maximum likelihood least squares identiﬁcation for systems with autoregressive moving average noise, Applied Mathematical Modelling 36 (5) (2012) 1842-1853. [34] J. Ding, F. Ding, X.P. Liu, G. Liu, Hierarchical least squares identiﬁcation for linear SISO systems with dual-rate sampled-data, IEEE Transactions on Automatic Control 56 (11) (2011) 2677-2683. [35] D.Q. Wang, F. Ding, Hierarchical least squares estimation algorithm for Hammerstein-Wiener systems, IEEE Signal Processing Letters 19 (12) (2012) 825-828. [36] D.Q. Wang, F. Ding, Least squares based and gradient based iterative identiﬁcation for Wiener nonlinear systems, Signal Processing 91 (5) (2011) 1182-1189. [37] F. Ding, Y.J. Liu, B. Bao, Gradient based and least squares based iterative estimation algorithms for multi-input multi-output systems, Proceedings of the Institution of Mechanical Engineers, Part I: Journal of Systems and Control Engineering 226 (1) (2012) 43-55. [38] F. Ding, X.P. Liu, G. Liu, Identiﬁcation methods for Hammerstein nonlinear systems, Digital Signal Processing 21 (2) (2011) 215-238. [39] D.Q. Wang, F. Ding, Y.Y. Chu, Data ﬁltering based recursive least squares algorithm for Hammerstein systems using the key-term separation principle, Information Sciences 222 (2013) 203-212. [40] F. Ding, G. Liu, X.P. Liu, Partially coupled stochastic gradient identiﬁcation methods for non-uniformly sampled systems, IEEE Transactions on Automatic Control 55 (8) (2010) 1976-1981. [41] D.Q. Wang, Least squares-based recursive and iterative estimation for output error moving average systems using data ﬁltering, IET Control Theory and Applications 5 (14) (2011) 1648-1657. [42] F. Ding, G. Liu, X.P. Liu, Parameter estimation with scarce measurements, Automatica 47 (8) (2011) 1646-1655. [43] F. Ding, Several multi-innovation identiﬁcation methods, Digital Signal Processing 20 (4) (2010) 1027-1039.

14