J. Appl. Math. Comput. https://doi.org/10.1007/s12190-017-1156-6 ORIGINAL RESEARCH

An iterative algorithm for solving the generalized Sylvester-conjugate matrix equation Caiqin Song1

Received: 24 August 2017 © Korean Society for Computational and Applied Mathematics 2017

Abstract This paper aims to extend the conjugate gradient least squares method to solve the least squares problem of the generalized Sylvester-conjugate matrix equation. For any initial values, the proposed iterative method can obtain the least squares Frobenius norm solution within finite iteration steps in the absence of roundoff errors. Finally, two numerical examples are given to show the efficiency of the presented iterative method. And it’s also proved that our proposed iterative algorithm is better than the existing LSQR iterative algorithm. Keywords Generalized Sylvester-conjugate matrix equations · Real inner · Least Frobenius norm solution Mathematics Subject Classification 65F10 · 65H10 · 15A33

1 Introduction Before starting this section, we first introduce the following notations which will be used in the rest of this paper. C m×n denotes the set of m × n complex matrices. For a matrix A ∈ C m×n , we denote its transpose, conjugate, conjugate transpose, trace, Frobenius norm and column space by A T , A, A H , tr(A), A and R(A), respectively. Let In and Sn denote the n×n unit matrix and reverse unit matrix respectively. The symbol vec(·) stands for the vec operator, i.e., for A = (a1 , a2 , . . . , an ) ∈ C m×n , where ai (i = 1, 2 . . . , n) denotes the ith column of A, vec(A) = (a1T , a2T , . . . , anT )T . For X

This project is granted financial support from NSFC (No. 11501246).

B 1

Caiqin Song [email protected] School of Mathematical Sciences, University of Jinan, Jinan 250022, China

123

C. Song

and Y two matrices in C m×n , we define real inner product as X, Y  = Re[tr(Y H X )]. The associated norm is the well-known Frobenius norm. For matrices R, A, B and X with appropriate dimension, a well-known property of the inner product is R, AX B = A H R B H , X . Two matrices X and Y are said to be orthogonal if X, Y  = 0. Many problems in mathematics, physics and engineering lead to solving linear matrix equations. For instance, a digital filter can be characterized by state variable equations x(n + 1) = Ax(n) + bu(n),

(1.1)

y(n) = cx(n) + du(n).

(1.2)

and

The solutions K and W of the Lyapunov matrix equations K = AK A T + bb T ,

(1.3)

W = AT W A + cT c

(1.4)

and

can analyzed the quantization noise generated by a digital filter [1]. The Lyapunov and (generalized) Sylvester matrix equations appear in robust control [2], neural network [3] and singular system control [4]. It can be shown that certain control problem, such as pole/eigenstructure assignment and observe design [5,6] of (1.1), are closely related with the following generalized Sylvester matrix equation p 

Ak X Bk +

k=1

q 

C j Y D j = E.

(1.5)

j=1

In this paper, we study the generalized Sylvester-conjugate matrix equation p  i=1

Ai X Bi +

t 

C j X D j = E,

(1.6)

j=1

with known matrices Ai ∈ C p×m , Bi ∈ C n×q , C j ∈ C p×m , D j ∈ C n×q , E ∈ C p×q and unknown matrix X ∈ C m×n for i = 1, 2, . . . , p, j = 1, 2, . . . , t. It is obvious that this class of matrix equation includes various linear matrix equations such as Lyapunov, Sylvester and the generalized Sylvester matrix equations as special cases. So far many techniques have been implemented for obtaining solutions of various linear matrix equations [7–17]. In [7], the least squares solution of the generalized Sylvester matrix equation AX B + CY D = E,

123

(1.7)

An iterative algorithm for solving the generalized...

was studied for symmetric arrowhead matrices with the least norm. By applying the hierarchical identification principle, Ding and Chen introduced gradient based iterative algorithms to solve (coupled) generalized Sylvester matrix equations, nonlinear systems [18–23]. Recently in [24–27], some efficient extended conjugate gradient methods were proposed for solving the least squares solution of several linear matrix equations. Very recently the matrix forms of CGS, BiCOR, CORS, Bi-CGSTAB and finite iterative methods were investigated for solving various linear matrix equations [28–37]. Based on the idea of the paper [7–9,22,26–37], this paper focuses on the goal of generalizing the CGLS algorithm for computing the least square solution of generalized Sylvester-conjugate matrix equation (1.6). We mainly consider the following two problems. Problem 1 For given Ai ∈ C p×m , Bi ∈ C n×q , C j ∈ C p×m , D j ∈ C n×q , E ∈ C p×q , i = 1, 2, . . . , p and j = 1, 2, . . . , t, find X ∗ ∈ C m×n such that     t   p  ∗ ∗ D − E  A X B + C X i i j j    i=1  j=1     p t    . = min  A X B + C X D − E i i j j   X ∈C m×n   i=1 j=1

(1.8)

Problem 2 Let G r denote the solution set of Problem 1, for given X 0 ∈ C m×n , find the matrix Xˇ ∈ Sr , such that  Xˇ − X 0  = min X − X 0 . X ∈G r

(1.9)

The rest of the paper is organized as follows. In Sect. 2, based on the CGLS method for solving Problems 1 and 2 we propose an iterative algorithm. It is shown that the proposed algorithm can obtain the solution of Problem 1 for any initial matrix within a finite number of iterations in the absence of roundoff errors. In Sect. 3, some numerical examples are given to show the efficiency of the proposed iterative method. An application example has been put in Sect. 4. Finally, we put some conclusions in Sect. 5. First of all, we introduce three lemmas which will be used in the next. And the proof line of Lemma 2 is similar to the reference [31]. Lemma 1 (Wang [38]) Let U be an inner product space, V be a subspace of U , and V ⊥ be the orthogonal complement subspace of V . For a given u ∈ U , if there exists a v0 ∈ V such that u − v0  ≤ u − v holds for any v ∈ V , then v0 is unique and v0 ∈ V is the unique minimization vector in V if and only if (u − v0 )⊥V, i.e. (u − v0 ) ∈ V ⊥ . Lemma 2 Assume that Rˆ is the residual of Eq. (1.6) corresponding to the matrix p  Xˆ ∈ C m×n , that is, Rˆ = E − i=1 Ai Xˆ Bi − tj=1 C j Xˆ D j . Then the matrix Xˆ is

123

C. Song

the least squares solution of Eq. (1.6) if p 

AiH Rˆ BiH +

i=1

t 

H H C j Rˆ D j = 0.

j=1

p Proof We first define the linear subspace W = {F|F = i=1 Ai X Bi + t p×m n×q p×m n×q , Bi ∈ C ,Cj ∈ C , Dj ∈ C , F ∈ C p×q , j=1 C j X D j } with Ai ∈ C p m×n and X ∈ C , i = 1, 2, . . . , p, j = 1, 2, . . . , t. Now if we let Fˆ = i=1 Ai Xˆ Bi + t ˆ ˆ j=1 C j X D j , then F ∈ W . Therefore, for any F ∈ W , we have  ˆ F = E − E − F,  ˆ = R,  =  =

p 

Ai Xˆ Bi −

t 

i=1 p 

t 

i=1

j=1

Ai X Bi +

p 



AiH Rˆ BiH ,

X +

i=1 p 

C j Xˆ D j ,

j=1



AiH Rˆ BiH +

i=1

Ai X Bi +

i=1

t 

 Cj X Dj

j=1

Cj X Dj  t 

 C j Rˆ D j , X H

j=1 t 

p 

H



C j Rˆ D j , X . H

H

j=1

p  H H Hence, if we let i=1 AiH Rˆ BiH + tj=1 C j Rˆ D j = 0, the above equations show ˆ F = 0. It follows from Lemma 1 that(E − F) ˆ ∈ W ⊥ . So the matrix Xˆ that E − F, is the least squares solution of Eq. (1.6).

Lemma 3 If we let Xˆ is a solution of Problem 1, then any solution X˜ of Problem 1 can be expressed as Xˆ + Z where the matrix Z ∈ C m×n satisfies p  i=1

Ai Z Bi +

t 

C j Z D j = 0.

(1.10)

j=1

Proof For any solution X˜ of Problem 1, if we define the matrix Z = X˜ − Xˆ , then we have X˜ = Xˆ + Z . Now it is showed that Eq. (1.10) holds. By applying Lemma 2, we can obtain 2    t    p   ˆ ˆ A + C − E X B X D i i j j     i=1 j=1  2   t   p   ˜ ˜ = Ai X Bi + C j X D j − E   i=1  j=1

123

An iterative algorithm for solving the generalized...

 2   t   p   ˆ ˆ = A ( X + Z )B + C ( X + Z )D − E i i j j    i=1  j=1  2   t   p   ˆ = Ai Z Bi + C j Z D j − R   i=1  j=1    p    t t     p  ˆ 2 = Ai Z Bi + Cj Z Dj Ai Z Bi + C j Z D j , Rˆ   +  R − 2  i=1  j=1 i=1 j=1      p  p t t      H H ˆ 2 = Ai Z Bi + Cj Z Dj AiH Rˆ BiH + C j Rˆ D j   + R − 2 Z ,  i=1  j=1 i=1 j=1    p  t    ˆ 2 = Ai Z Bi + Cj Z Dj   +  R .  i=1  j=1 This shows that the Eq. (1.10) holds.



2 Main results The CGLS method is a powerful method for solving the solution of large sparse least squares problem minn b − Ax. In this section, we first propose an iterative X ∈C

algorithm by generalizing the CGLS method to solve Problem 1, then present some basic properties of the algorithm. We also consider finding the least Frobenius norm solution of Problem 1. Algorithm CGLS. [39] x(1) ∈ R n is an initial guess, r (1) = b − Ax(1), s(1) = A T r (1), p(1) = s(1), γ (1) = s(1)2 , for k = 1, 2, 3, . . . repeat the following: q(k) = Ap(k) δ(k) = γ (k)/q(k)2 x(k + 1) = x(k) + δ(k) p(k) r (k + 1) = r (k) − δ(k)q(k) s(k + 1) = A T r (k + 1) γ (k + 1) = s(k + 1)2 λ(k) = γ (k + 1)/γ (k) p(k + 1) = s(k + 1) + λ(k) p(k). Based on the above CGLS method, we propose the following matrix form of CGLS iterative algorithm (CGLS-M) to solve the least squares problem of the generalized Sylvester-conjugate matrix equation (1.6).

123

C. Song

Algorithm 1 Step 1. Input matrices Ai ∈ C p×m , Bi ∈ C n×q , C j ∈ C p×m , D j ∈ C n×q , E ∈ C p×q and X (1) ∈ C m×n for i = 1, 2, . . . , p, j = 1, 2, . . . , t; Step 2. Compute R(1) = E −

p 

Ai X (1)Bi −

i=1

S(1) =

p 

t 

C j X (1)D j ;

j=1

AiH R(1)BiH +

i=1

t 

H

H

C j R(1)D j ;

j=1

P(1) = S(1), γ (1) = S(1)2 ; For k = 1, 2, 3, . . . repeat the following: Step 3. If R(k) = 0, then stop and X (k) is the solution of Eq. (1.6), break; else if R(k) = 0 but S(k) = 0, then stop and X (k) is the solution of Problem 1, break; else k := k + 1; Step 4. Compute Q(k) =

p 

Ai P(k)Bi +

i=1

t 

C j P(k)D j ;

j=1

δ(k) = γ (k)/Q(k)2 ; X (k + 1) = X (k) + δ(k)P(k); R(k + 1) = R(k) − δ(k)Q(k); t p H H AiH R(k + 1)BiH + C j R(k + 1)D j S(k + 1) = i=1 j=1   p t H H ; = S(k) − δ(k) AiH Q(k)BiH + C j Q(k)D j i=1

j=1

γ (k + 1) = S(k + 1)2 ; λ(k) = γ (k + 1)/γ (k); P(k + 1) = S(k + 1) + λ(k)P(k). Step 5. Go to Step 3. Some basic properties of Algorithm 1 are listed in the following lemmas. Lemma 4 For the sequences S(k), P(k) and Q(k) which are generated by Algorithm 1, if there exists a positive number r , such that S(u) = 0 and Q(u) = 0, ∀u = 1, 2, . . . , r , then the following statements hold for u, w = 1, 2, . . . , r and u = w. (1) S(u), S(w) = 0, (2) Q(u), Q(w) = 0, (3) P(u), S(w) = 0.

123

An iterative algorithm for solving the generalized...

Proof Step 1: We prove the conclusion by induction. Because the real inner product is commutative, it is enough to prove three statements for 1 ≤ u < w ≤ r. For u = 1, w = 2, by Algorithm 1, we have  p 

t H H S(1), S(2) = S(1), S(1) − δ(1) AiH Q(1)BiH + C j Q(1)D j i=1 j=1

p t H H 2 H H = S(1) −δ(1) S(1), Ai Q(1)Bi + C j Q(1)D j i=1 j=1  p

t 2 Ai S(1)Bi + C j S(1)D j , Q(1) = S(1) − δ(1) i=1

j=1

= S(1) − δ(1)Q(1), Q(1) = 0, (2.1)

p t Q(1), Q(2) = Q(1), Ai P(2)Bi + C j P(2)D j ) i=1 j=1 p Ai (S(2) + λ(1)P(1))Bi = Q(1), i=1

t + C j (S(2) + λ(1)P(1))D j j=1

p t = λ(1)Q(1)2 + Q(1), Ai S(2)Bi + C j S(2)D j i=1 j=1  p 1 A H (R(2) − R(1))BiH = λ(1)Q(1)2 − i=1 i δ(1)

t H H + C j (R(2) − R(1))D j , S(2) 2

j=1

1 S(2) − S(1), S(2) = 0, = λ(1)Q(1)2 − δ(1)

(2.2)

and P(1), S(2) = S(1), S(2) = 0.

(2.3)

Step 2: In this case, for u < w < r we assume that S(u), S(w) = 0,

Q(u), Q(w) = 0,

P(u), S(w) = 0.

Thus, we can write S(u), S(w + 1)  p 

t H H = S(u), S(w) − δ(w) AiH Q(w)BiH + C j Q(w)D j i=1 j=1

p t H H = −δ(w) S(u), AiH Q(w)BiH + C j Q(w)D j i=1 j=1  p

t = −δ(w) Ai S(u)Bi + C j S(u)D j , Q(w) i=1 j=1  p Ai (P(u) − λ(u − 1)P(u − 1))Bi = −δ(w) i=1

123

C. Song

t +

j=1

C j (P(u) − λ(u − 1)P(u − 1))D j , Q(w)



= −δ(w)Q(u) − λ(u − 1)Q(u − 1), Q(w) = 0, (2.4) Q(u), Q(w + 1)

t p Ai P(w + 1)Bi + C j P(w + 1)D j = Q(u), i=1 j=1 p Ai (S(w + 1) + λ(w)P(w))Bi = Q(u), i=1

t + C j (S(w + 1) + λ(w)P(w))D j j=1

p t = Q(u), Ai S(w + 1)Bi + C j S(w + 1)D j + λ(w)Q(w) i=1 j=1

p t 1 R(u + 1) − R(u), Ai S(w + 1)Bi + C j S(w + 1)D j =− i=1 j=1 δ(u) 1  p =− A H (R(u + 1) − R(u))BiH i=1 i δ(u)

t H H + C j (R(u + 1) − R(u))D j , S(w + 1) j=1

1 =− S(u + 1) − S(u), S(w + 1) = 0, δ(u) P(u), S(w + 1)

 p t H H AiH Q(w)BiH + C j Q(w)D j = P(u), S(w) − δ(w) i=1 j=1

p t H H H H = −δ(w) P(u), Ai Q(w)Bi + C j Q(w)D j i=1 j=1  p

t = −δ(w) Ai P(u)Bi + C j P(u)D j , Q(w) = 0. i=1

j=1

(2.5)

(2.6)

For u = w, we have S(w), S(w + 1)  p 

t H H = S(w), S(w) − δ(w) AiH Q(w)BiH + C j Q(w)D j i=1 j=1

p t H H = S(w)2 − δ(w) S(w), AiH Q(w)BiH + C j Q(w)D j i=1 j=1  p

t 2 = S(w) − δ(w) Ai S(w)Bi + C j S(w)D j , Q(w) i=1 j=1  p 2 Ai (P(w) − λ(w − 1)P(w − 1))Bi = S(w) − δ(w) i=1

t + C j (P(w) − λ(w − 1)P(w − 1))D j , Q(w) j=1

= S(w) − δ(w)Q(w) − λ(w − 1)Q(w − 1), Q(w) = 0, Q(w), Q(w + 1)

p t = Q(w), Ai P(w + 1)Bi + C j P(w + 1)D j 2

i=1

123

j=1

(2.7)

An iterative algorithm for solving the generalized...

p = Q(w), Ai (S(w + 1) + λ(w)P(w))Bi i=1

t + C j (S(w + 1) + λ(w)P(w))D j j=1

p t = Q(w), Ai S(w + 1)Bi + C j S(w + 1)D j ) + λ(w)Q(w) i=1 2

j=1

= λ(w)Q(w)

p t 1 R(w + 1) − R(w), − Ai S(w + 1)Bi + C j S(w + 1)D j i=1 j=1 δ(w)  p 1 = λ(w)Q(w)2 − A H (R(w + 1) − R(w))BiH i=1 i δ(w)

t H H + C j (R(w + 1) − R(w))D j , S(w + 1) j=1

1 S(w + 1) − S(w), S(w + 1) = λ(w)Q(w)2 − δ(w) 1 = λ(w)Q(w)2 − S(w + 1)2 = 0, δ(w) P(w), S(w + 1)

 p t H H AiH Q(w)BiH + C j Q(w)D j = P(w), S(w) − δ(w) i=1

(2.8)

j=1

= S(w) + λ(w − 1)S(w − 1) + λ(w − 1)λ(w − 2)S(w − 2) + · · · + λ(w − 1) . . . λ(1)S(1), S(w)

t p H H AiH Q(w)BiH + C j Q(w)D j − δ(w) P(w), i=1 j=1  p

t 2 Ai P(w)Bi + C j P(w)D j , Q(w) = 0. (2.9) = S(w) − δ(w) i=1

j=1



By the principle of induction, we draw the conclusion.

Lemma 5 If there exists a positive number l such that δ(l) = 0 or δ(l) = ∞ in Algorithm 1, then X (l) is a solution of Problem 1. Proof If δ(l) = 0, we have S(l)2 = 0. If δ(l) = ∞, we have Q(l)2 = 0. It is showed that S(l)2 = S(l) + λ(l − 1)S(l − 1) + λ(l − 1)λ(l − 2)S(l − 2) + · · · + λ(l − 1) . . . λ(1)S(1), S(l) = P(l), S(l) p t = P(l), AiH R(l)BiH + i=1

j=1

H

C j R(l)D j

H



= 0. Therefore, for δ(l) = 0 or δ(l) = ∞ we have S(l) =

p  i=1

AiH R(l)BiH +

t 

H

C j R(l)D j

H

= 0.

j=1

123

C. Song

So by Lemma 2 we can conclude that X (l) is the solution of Problem 1.



Theorem 1 If the matrix equation (1.6) is consistent, then, for any arbitrary initial matrix X (l), the solution X ∗ of the Problem 1 can be obtained by using Algorithm 1 within a finite number of iterations in the absence of roundoff errors. Proof By Lemma 4, the set S(i), i = 1, 2, 3, . . . , m × n is an orthogonal basis of the real inner product space C m×n with dimension m × n. Therefore, we can obtain S(m × n + 1) = 0. It is also showed that X (m × n + 1) is the solution of Problem 1 in the absence of roundoff errors.

p t H H H Theorem 2 Supposed that the initial matrix is j=1 C j i=1 Ai F(l)Bi + H

F(l)D j , where F(1) ∈ C m×n is an arbitrary matrix, or especially X (1) = 0, then the solution X ∗ generated by Algorithm 1 is the least norm solution of Problem 1. p  H H Proof Let the initial matrix be X (1) = i=1 AiH F(l)BiH + tj=1 C j F(l)D j , then, it follows from Algorithm 1 that the generated matrix X (k) can be easily  p H H expressed as X (k) = i=1 AiH F(k)BiH + tj=1 C j F(k)D j for certain matrices F(k) ∈ C m×n for k = 2, 3, . . .. This implies that there exists a matrix F ∗ ∈ C m×n p  H ∗ H such that X ∗ = i=1 AiH F ∗ BiH + tj=1 C j F D j . Now if we let X˜ ∗ is an arbitrary solution of Problem 1, then, by using Lemma 3, there exists the matrix Z ∗ ∈ C m×n p ∗ ∗ ∗ ∗ such that X˜ = X + Z and i=1 Ai Z Bi + tj=1 C j Z ∗ D j = 0. Thus, one can obtain  p  t   H ∗ H ∗ ∗ H ∗ H ∗ X , Z  = Ai F Bi + Cj F Dj , Z 

i=1 ∗

= F ,

j=1 p 

∗

Ai Z Bi +

i=1

t 

 Cj

Z∗D

j

= 0.

(2.10)

j=1

By applying Eq. (2.10), one can obtain  X˜ ∗ 2 = X ∗ + Z ∗ 2 = X ∗ 2 + Z ∗ 2 + 2X ∗ , Z ∗  = X ∗ 2 + Z ∗ 2 ≥ X ∗ 2 This shows that the solution X ∗ is the least Frobenius norm solution of Problem 1.

Similar to [31], the minimization property of the proposed algorithm is stated as follows. This property shows that Algorithm 1 converges smoothly. Theorem 3 For any initial matrix X (1) ∈ C m×n , we have 2    t    p   A X (k + 1)B + C X (k + 1)D − E i i j j     i=1 j=1 2    t    p  = min  Ai X Bi + C j X D j − E  , X ∈ψk   i=1 j=1

123

(2.11)

An iterative algorithm for solving the generalized...

where X (k + 1) is generated by Algorithm 1 at the k + 1-th iteration and ψk presents an affine subspace which has the following form ψk = X (1) + span(P(1), P(2), . . . , P(k)).

(2.12)

Proof For any matrix X ∈ ψk , it follows from Eq. (2.12) that there exist numbers α1 , α2 , . . . , αk such that X = X (1) +

k 

αl P(l).

(2.13)

l=1

Now we define the continuous and differentiable function f with respect to the variable α1 , α2 , . . . , αk as f (α1 , α2 , . . . , αk )  2  s  k t k       = Ai (X (1) + αl P(l))Bi + C j (X (1) + αl P(l))D j − E   .  i=1  l=1 j=1 l=1 By Lemma 4 one can obtain s  f (α1 , α2 , . . . , αk ) = 

i=1

Ai X (1)Bi +

s k αl + l=1

= R(1)2 +

k

i=1

t j=1

C j X (1)D j − E

Ai P(l)Bi +

α 2 Q(1)2 l=1 l

 2 C j P(l)D j   j=1

t

− 2αl Q(1), R(1).

Now we consider the problem of minimizing the function f (α1 , α2 , . . . , αk ). It is obvious that 2    p t     min f (α1 , α2 , . . . , αk ) = min  Ai X Bi + C j X D j − E  . αl X ∈ψk   i=1 j=1 For this function, the minimum occurs when ∂(α1 , α2 , . . . , αk ) = 0 for l = 1, 2, . . . , k. ∂αl Thus, we can get αl =

Q(l), R(l) . Q(l)2

123

C. Song

From Algorithm 1, one can obtain R(1) = R(l) + δ(l − 1)Q(l − 1) + δ(l − 2)Q(l − 2) + · · · + δ(1)Q(1). Therefore, it follows from Lemma 4 that

s  H H P(l), i=1 AiH R(l)BiH + tj=1 C j R(l)D j Q(l), R(l) αl = = Q(l)2 Q(l)2 P(l), S(l) S(l) + λ(l − 1)P(l − 1), S(l) = = 2 Q(l) Q(l)2 2 S(l) = = δ(l). Q(l)2



Thus, we complete the proof.

By Theorem 1, the solution generalized by Algorithm 1 at the k + 1−th iteration for any initial matrix minimizes the residual norm in the affine subspace ψk . Also one has 2    t    s  Ai X (k + l)Bi + C j X (k + l)D j − E      i=1 j=1  2   t   s   , ≤ A X (k)B + C X (k)D − E i i j j    i=1  j=1

(2.14)

which shows that the sequence of the norm of residuals R(1), R(2), . . . is monotonically decreasing. This decent property of the norm of the residuals shows that the Algorithm 1 processes fast and smoothly convergence. Now we will solve Problem 2. For a given matrix X 0 , by Problem 1 we can get     t    p  min  A X B + C X D − E i i j j   m×n X ∈C   i=1 j=1   p t   = min  A (X − X )B + C j (X − X 0 )D j i 0 i  X ∈C m×n  i=1 j=1 ⎞ ⎛  p t    − ⎝E − Ai X 0 Bi − C j X 0 D j ⎠ .  i=1 j=1 If we let the set E 1 = E −

p  i=1

Ai X 0 Bi −

t  j=1

C j X 0 D j and X 1 = X − X 0 , then

Problem 2 is equivalent to find the least Frobenius norm solution X 1∗ of

123

(2.15)

An iterative algorithm for solving the generalized...

    t   p   min  Ai X 1 Bi + C j X 1 D j − E1 , X 1 ∈C m×n   i=1 j=1

(2.16)

which can be computed by using Algorithm 1 with the initial matrix X 1 (1) = s t H H H H m×n is an arbitrary matrix, or espei=1 Ai F Bi + j=1 C j F D j where F ∈ C cially X 1 (1) = 0. Thus the solution of Problem 2  Xˇ − X 0  = min X − X 0  can X ∈G r

be stated as Xˇ = X 1∗ + X 0 .

(2.17)

It’s showed that based on the Problem 1, we can solve the Problem 2. In other words, we can obtain the solution of Problem 2 by using the solution of Problem 1.

3 Numerical examples In this section, we shall give some numerical examples to illustrate the efficiency of Algorithm 1. All the tests are performed by MATLAB with machine precision around 10−16 . Example 3.1 Find the least Frobenius norm solution of the following matrix equation AX B + C X D = M,

(3.1)

where ⎛

⎞ ⎛ ⎞ 1 + i 1 − i 6i 4 + 3i 2 3 − 2i 4 + 6i 9 + 8i ⎜2 + i 0 ⎟ ⎜ ⎟ 12 10 12 9 ⎟ , B = ⎜ 4 + 6i 11 ⎟, A=⎜ ⎝ 5 + 6i 2 + 3i 11 − 2i i ⎠ ⎝0 ⎠ 12 15 18 1 12 0 9i 2 − i −9i 12 11 ⎞ ⎛ ⎛ ⎞ 9 + 2i −i i 6 + 8i i 9 15 2i ⎜ ⎜ 11 19 0 11 ⎟ 2 + 3i 12 11 ⎟ ⎟ , D = ⎜ 4i ⎟ C=⎜ ⎝ 23 26 + i 9i 9 ⎠ , ⎝ 18 11 − 9i 21 i ⎠ 9 + 8i 11 i 8i 23 0 6 8i ⎞ ⎛ 4660 − 5011i 6322 − 1536i 2161 + 1789i 6035 + 1940i ⎜ 14900 − 1787i 13144 − 1058i 5872 + 5983i 12675 + 8798i ⎟ ⎟ M=⎜ ⎝ 12135 − 12509i 12910 − 5709i 4671 + 6089i 11126 + 4083i ⎠ . 10031 + 8092i 12415 + 3909i 3793 + 9126i 5328 + 12279i It can be verified that the generalized Sylvester-conjugate matrix equation (3.1) are consistent and have the solution ⎛ ⎞ 4 + 3i 2 + i 11 6 ⎜ 9 + 6i 11 + 2i 8 + 9i 9 ⎟ ⎟. X =⎜ ⎝ 9i 4 2i 11i ⎠ 23 4 7 12

123

C. Song 5

0

r(k)

−5

−10

−15

CGLS− M LSQR−M

−20

0

20

40

60

80

100

120

140

160

k(iteration number) Fig. 1 The residual of solution for Example 3.1

If we let the initial matrix X 1 = 0, applying Algorithm 1, we obtain the solution ⎛

4.0000 + 2.9999i ⎜ 8.9999 + 6.0000i X (126) = ⎜ ⎝ 8.9999i 22.9999

1.9999 + 0.99999i 11.0000 + 2.0000i 3.9999 4.0000

11.0000 7.9999 + 9.0000i 1.9999i 6.9999

⎞ 5.9999 9.0000 ⎟ ⎟. 11.0000i ⎠ 12.0000

with corresponding residual r (126) = 9.8915 × 10−16 , and relative error e(126) = 1.3293 × 10−14 . The obtained results are presented in Figs. 1 and 2, where rk = log10 E − AX (k)B − C X (k)D, and err (k) =

X (k) − X  . X 

If we let the initial matrix X (1) = 0, X (1) = I4 , X (1) = ones(4, 4), respectively, applying Algorithm 1, the residual and the relative error of solution are presented in Figs. 3 and 4.

123

An iterative algorithm for solving the generalized... 0 -2 LSQR-M

e(k)

-4

CGLS-M

-6 -8

-10 -12 -14

0

20

40

60

80

100

120

140

160

k(iterative number)

Fig. 2 The relative error of solution for Example 3.1 5

0

r(k)

-5

-10 X 1 =I

4

X 1 =ones(4,4)

-15

X 1 =zeros(4,4)

-20

0

20

40

60

80

100

120

140

k(iterative number) Fig. 3 The residual of solution for Example 3.1

Remark 1 From Figs. 1 and 2, the residual and the relative error of the solution to matrix equation (3.1) are reduced with the addition of iterative number. It’s showed that the iterative solution converges to the exact solution. Moreover, with the addition of iterative number, the residual and the relative error of LSQR-M iterative algorithm are not variety. Finally, the residual and the relative error of iterative algorithm are smaller than those of LSQR-M iterative algorithm. It is proved that our proposed method (CGLS-M iterative algorithm) is efficient and is a good algorithm.

123

C. Song 0

X =zeros(4,4) 1

X =ones(4,4) 1

X =I

-5

4

e(k)

1

-10

-15

0

20

40

60

80

100

120

140

k(iterative number) Fig. 4 The relative error of solution for Example 3.1

Remark 2 From Figs. 3 and 4, for the different initial value X (1) = 0, X (1) = I4 , X (1) = ones(4, 4), respectively, the iterative solution to the matrix equation (3.1) all converges to the exact solution to the matrix equation (3.1). Example 3.2 Find the least Frobenius norm solution of the following matrix equation AX + X D = M,

(3.2)

where ⎞ 1 + 2i 11 − 3i 12 + 6i 4 + 23i 11 2 ⎟ ⎜ 2 + 11i 0 12 10 34 + 5i 11 ⎟ ⎜ ⎟ ⎜ 5 + 6i 2 + 3i 11 − 2i i 34 45 ⎟, ⎜ A=⎜ ⎟ 1 12 0 9i 13 + 45i 3 ⎟ ⎜ ⎝ 23 43 21 99i 45i 4 + 9i ⎠ 34 89 45i 2 + 3i 4 + 7i 4 ⎛ ⎞ 9 + 2i −i i 6 + 8i 12 23 ⎜ 4i 19 0 11 3i 3 + 4i ⎟ ⎜ ⎟ ⎜ 23 ⎟ 26 + i 9i 9 12 23 ⎜ ⎟, D=⎜ 0 6 8i 2 18 + 9i ⎟ ⎜ 23 ⎟ ⎝ 112 123 + i 119 21 ⎠ 29 + 9i 3i 11 2 3i 4 + 28i 89 + 9i 788 ⎛

⎛

1349 + 1310i ⎜ 2951 + 2863i ⎜ ⎜ 574 + 3690i M =⎜ ⎜ −802 + 1004i ⎜ ⎝ 14159 + 2417i 24695 − 1540i

123

829 − 172i 4222 + 22i 4159 − 959i 1885 − 105i 16561 + 1277i 29322 + 289i

604 + 229i 4948 + 311i 6374 − 1083i 2110 + 827i 16974 + 3140i 25160 + 2003i

1547 + 2029i 5997 + 782i 5778 + 166i 2515 + 7660i 2939 + 6494i 14052 − 275i

8166 + 846i 7743 − 2175i 14398 − 1699i 10150 + 6537i 6763 + 9222i 12010 − 18641i

⎞ 53917 + 3054i 2060 − 23329i ⎟ ⎟ 10406 − 3223i ⎟ ⎟. 62506 + 2197i ⎟ ⎟ 11737 + 10689i ⎠ 4965 − 182152i

An iterative algorithm for solving the generalized... 5

0

r(k)

-5

-10

CGLS-M

-15

-20

LSQR-M

0

100

200

300

400

500

600

k( iteration number )

Fig. 5 The residual of solution for Example 3.2

It can be verified that the generalized Sylvester-conjugate matrix equation (3.2) are consistent and have the solution ⎛

4 + 3i ⎜ 9 + 6i ⎜ ⎜ 9i X =⎜ ⎜ 23 ⎜ ⎝ 2 + 78i 56i

2+i 11 + 2i 4 4 2 + 9i 79

11 8 + 9i 2i 7 17 118

6 9 11i 12i 119 19

3i 23 3 + 10i 12 + 3i 127 191

⎞ 67 ⎟ 34i ⎟ 12 + 18i ⎟ ⎟. ⎟ 78 ⎟ ⎠ 12 237i

If we let the initial matrix X (1) = 0, applying Algorithm 1, we obtain the solutions, that are ⎛

4.0000 + 3.0000i ⎜ 8.9999 + 6.0000i ⎜ ⎜ 8.9999i X (587) = ⎜ ⎜ 22.9999 ⎜ ⎝ 2.0000 + 77.9999i 56.0000i 6.0000 9.0000 11.0000i 12.0000i 118.9999 19.0000

2.0000 + 0.99999i 10.9999 + 2.0000i 3.9999 3.9999 2.0000 + 8.9999i 79.0000

2.9999i 23.0000 3.0000 + 1.0000i 12.0000 + 3.0000i 126.9999 191.0000

11.0000 7.9999 + 9.0000i 1.9999i 6.9999 17.0000 117.9999 ⎞

67.0000 ⎟ 33.9999i ⎟ 12.0000 + 18.0000i ⎟ ⎟. ⎟ 78.0000 ⎟ ⎠ 11.9999 237.0000

123

C. Song 0

CGLS-M

-5

e(k)

LSQR-M

-10

-15

0

100

200

300

400

500

600

500

600

k(iterative number) Fig. 6 The relative error of solution for Example 3.2 5

0

r(k)

-5

X =I

-10

1

6

X 1 =ones(6,6) X =zeros(6,6) 1

-15

-20

0

100

200

300

400

k(iterative number) Fig. 7 The residual of solution for Example 3.2

with corresponding residual r (587) = 6.9270 × 10−16 . and relative error e(587) = 7.9381 × 10−14 .

123

An iterative algorithm for solving the generalized... 0

X 1 =zero(6,6) X 1 =I 6

-5

e(k)

X 1 =ones(6,6)

-10

-15

0

100

200

300

400

500

600

k(iterative number) Fig. 8 The relative error of solution for Example 3.2

The obtained results are presented in Figs. 5 and 6,where rk = log10 E − AX (k)B − C X (k)D, and err (k) =

X (k) − X  . X 

If we let the initial matrix X (1) = 0, X (1) = I4 , X (1) = ones(4, 4), respectively, applying Algorithm 1, the residual and the relative error of solution are presented in Figs. 7 and 8. Remark 3 From Figs. 5 and 6, it’s showed that the iterative solution converges to the exact solution. Moreover, with the addition of iterative number, the residual and the relative error of iterative algorithm are smaller than those of LSQR-M iterative algorithm. It is proved that our proposed method (CGLS-M iterative algorithm) is efficient and is a good algorithm. From Figs. 7 and 8, for the different initial value X (1) = 0, X (1) = I4 , X (1) = ones(4, 4), respectively, the iterative solution to the matrix equation (3.2) all converges to its exact solution.

4 Application example In linear systems, the problem of eigenstructure assignment is to find a control law such that the system matrix of the closed-loop system has the desired eigenvalues and multiplicities, and simultaneously determine the corresponding eigenvectors and generalized eigenvectors.

123

C. Song

Now we consider the following discrete-time anti-linear system x(t + 1) = Ax(t) + Bu(t)

(4.1)

where x(t) ∈ C n is the state, and u(t) ∈ C r is the input. A ∈ C n×n and C ∈ C n×r are the system matrix and input matrix, respectively. If the state feedback u(t) = K x(t),

(4.2)

is applied to the system (4.1), the closed-loop system is obtained as x(t + 1) = (A + B K )x(t),

(4.3)

Given a prescribed matrix F ∈ C n×n , find a state feedback controller in the form of (4.2) such that the system matrix (A + B K ) of the closed-loop system (3) is consimilar to the matrix F, and simultaneously determine a nonsingular matrix X ∈ C n×n satisfying X −1 (A + B K )X = F,

(4.4)

Due to the non-singularity of X , the equation (4.4) can equivalently written as AX + B K X = X F,

(4.5)

Let K X = W , then, the equation (4.5) can be transformed into the following equation AX + BW = X F,

(4.6)

which is the so-called con-Sylvester matrix equation. By Algorithm 1, we can obtain the solution of equation (4.6).

5 Conclusions In this paper, the CGLS iterative algorithm is constructed to solve the least Frobenius norm solution of generalized Sylvester-conjugate matrix equation (1.6). By this algorithm, for any initial matrix X (1), a solution X ∗ can be obtained in finite iteration steps in the absence of roundoff errors, and the least Frobenius norm solution can be obtained by choosing any initial matrix. In addition, by use of this iterative method, the optimal approximation solution Xˇ to a given matrix X 0 can be derived by first finding the least Frobenius norm solution of a new corresponding matrix equation. And we have checked the obtained theory by the numerical examples. The given numerical examples shows that the proposed iterative algorithms are quite efficient than the existing LSQR iterative algorithm [25–27]. Acknowledgements The authors are very grateful to the anonymous reviewers and the editor for their helpful comments and suggestions which have helped us in improving the quality of this paper.

123

An iterative algorithm for solving the generalized...

References 1. Gerheim, A.: Numerical solution of the Lyapunov equation for narrow-band digital filters. IEEE Trans. Circuits Syst. 31(11), 991–992 (1984) 2. Varga, A.: Robust pole assignment via Sylvester equation based state feedback parametrization. In: Proceedings of the 2000 IEEE International Symposium On Computer-aided Control System Design, Alsaka, USA (2000) 3. Zhang, Y.N., Jiang, D.C., Wang, J.: A recurrent neural network for solving Sylvester equation with time-varying coefficients. IEEE Trans. Neural Netw. 13, 1053–1063 (2002) 4. Shahzad, A., Jones, B.L., Kerrigan, E.C., Constantinides, G.A.: An efficient algorithm for the solution of a coupled Sylvester equation appearing in descriptor systems. Automatica 47, 244–248 (2011) 5. Duan, G.R.: Parametric eigenstructure assignment in high-order linear systems. Int. J. Control Autom. Syst. 3, 419–429 (2005) 6. Zhou, B., Duan, G.R.: On the generalized Sylvester mapping and matrix equationa. Syst. Control Lett. 57, 200–208 (2008) 7. Li, F.L., Hu, X.Y., Zhang, L.: The generalized anti-reflexive solutions for a class of matrix equations B X = C, X D = E. Comput. Appl. Math. 27, 31–46 (2008) 8. Yuan, Y.X., Dai, H.: Generalized reflexive solutions of the matrix equation AX B = D and associated optimal approximation problem. Comput. Math. Appl. 56, 1643–1649 (2008) 9. Zhang, J.C., Zhou, S.Z., Hu, X.Y.: The (P, Q)generalized reflexive and anti-reflexive solutions of the matrix equation AX = B. Appl. Math. Comput. 209, 254–258 (2009) 10. Shen, X.P., Chen, G.L.: An iterative method for the symmetric and skew symmetric solutions of a linear matrix equation AX B + CY D = E. J. Comput. Appl. Math. 233, 3030–3040 (2010) 11. Wang, Q.W.: Bisymmetric and centrosymmetric solutions to systems of real quaternion matrix equations. Comput. Math. Appl. 49, 641–650 (2005) 12. Wang, Q.W.: The general solution to a system of real quaternion matrix equations. Comput. Math. Appl. 49, 665–675 (2005) 13. Wang, Q.W., Li, C.K.: Ranks and the least-norm of the general solution to a system of quaternion matrix equations. Linear Algebra Appl. 430, 1626–1640 (2009) 14. Wang, Q.W., He, Z.H.: Systems of coupled generalized Sylvester matrix equations. Automatica 50, 2840–2844 (2014) 15. Wang, Q.W., He, Z.H.: Solvability conditions and general solution for the mixed Sylvester equations. Automatica 49, 2713–2719 (2013) 16. Sheng, X.P., Chen, G.L.: A finite iterative method for solving a pair of linear matrix equations (AX B, C X D) = (E, F). Appl. Math. Comput. 189, 1350–1358 (2007) 17. Wang, M.H., Cheng, X.H., Wei, M.S.: Iterative algorithms for solving the matrix equation AX B + C X T D = E. Appl. Math. Comput. 187, 622–629 (2007) 18. Ding, F., Chen, T.: Iterative least squares solutions of coupled Sylvester matrix equations. Syst. Control Lett. 54, 95–107 (2005) 19. Ding, F., Chen, T.: On iterative solutions of general coupled matrix equations. SIAM J. Control. Optim. 44, 2269–2284 (2006) 20. Ding, F., Chen, T.: Gradient based iterative algorithms for solving a class of matrix equations. IEEE Trans. Autom. Control 50, 1216–1221 (2005) 21. Ding, F., Wang, Y., Ding, J.: Recursive least squares parameter identification algorithms for systems with colored noise using the filtering technique and the auxilary model. Digit. Signal Proc. 37, 100–108 (2015) 22. Ding, F., Wang, X., Chen, Q., Xiao, Y.: Recursive least squares parameter estimation for a class of output nonlinear systems based on the model decomposition. Signal Process Circuits Syst (2016). https://doi.org/10.1007/s00034-015-0190-6 23. Li, H., Gao, Z., Zhao, D.: Least squares solutions of the matrix equation AX B + CY D = E with the least norm for symmetric arrowhead matrices. Appl. Math. Comput. 226, 719–724 (2014) 24. Hajarian, M.: Developing the CGLS algorithm for the least squares solutions of the general coupled matrix equations. Math. Methods Appl. Sci. 37, 2782–2798 (2014) 25. Peng, Z.: The reflexive least squares solutions of the matrix equation A1 X 1 B1 + A2 X 2 B2 + · · · + Al X l Bl = C with a submatrix constraint. Numer. Algorithms 64, 455–480 (2013) 26. Peng, Z., Xin, H.: The reflexive least squares solutions of the general coupled matrix equations with a submatrix constraint. Appl. Math. Comput. 225, 425–445 (2013)

123

C. Song 27. Peng, Z.: The (R, S)-symmetric least squares solutions of the general coupled matrix equations. Linear Multilinear Algebra 63, 1086–1105 (2015) 28. Hajarian, M.: Developing BiCOR and CORS methods for coupled Sylvester-transpose and periodic Sylvester matrix equations. Appl. Math. Model. 39, 6073–6084 (2015) 29. Hajarian, M.: Matrix iterative methods for solving the Sylvester-transpose and periodic Sylvester matrix equations. J. Franklin Inst. 350, 3328–3341 (2013) 30. Hajarian, M.: Matrix form of the CGS method for solving general coupled matrix equations. Appl. Math. Lett. 34, 37–42 (2014) 31. Hajarian, M.: Extending the CGLS algorithm for least squares solutions of the generalized Sylvestertranspose matrix equations. J. Franklin Inst. 353, 11688–1185 (2016) 32. Dmytryshyn, A., Kastrom, B.: Coupled Sylvester-type matrix equations and block diagonalization. SIAM J. Matrix Anal. Appl. 38, 580–593 (2015) 33. Wu, A.G., Zhang, E.Z., Liu, F.C.: On closed-form solutions to the generalized Sylvester-conjugate matrix equation. Appl. Math. Comput. 218(19), 9730–9741 (2012) 34. Wu, A.G., Lv, L.L., Duan, G.R., Liu, W.Q.: Parametric solutions to Sylvester-conjugate matrix equations. Comput. Math. Appl. 62(12), 4806–4806 (2011) 35. Wu, A.G., Lv, L.L., Hou, M.Z.: Finite iterative algorithms for extended Sylvester-conjugate matrix equations. Math. Comput. Modell. 54(9–10), 2363–2384 (2011) 36. Wu, A.G., Lv, L.L., Li, B., Zhang, Y., Duan, G.R.: Finite iterative solutions to coupled Sylvesterconjugate matrix equations. Appl. Math. Model. 35(3), 1065–1080 (2011) 37. Wu, A.G., Sun, Y., Feng, G.: Closed-form solution to the non-homogeneous generalised Sylvester matrix equation. IET Control Theory Appl. 4(10), 1914–1921 (2010) 38. Wang, R.C.: Functional Analysis and Optimization Theory. Beijing University of Aeronautics and Astronautics Press, Beijing (2003) 39. Björck, Å.: Numerical Methods for Least Squares Problems. SIAM, Philadelphia (1996)

123