Processing math: 84%
Research article

Conjugate gradient algorithm for consistent generalized Sylvester-transpose matrix equations

  • Received: 26 October 2021 Revised: 13 December 2021 Accepted: 21 December 2021 Published: 06 January 2022
  • MSC : 15A60, 15A69, 65F45

  • We develop an effective algorithm to find a well-approximate solution of a generalized Sylvester-transpose matrix equation where all coefficient matrices and an unknown matrix are rectangular. The algorithm aims to construct a finite sequence of approximated solutions from any given initial matrix. It turns out that the associated residual matrices are orthogonal, and thus, the desire solution comes out in the final step with a satisfactory error. We provide numerical experiments to show the capability and performance of the algorithm.

    Citation: Kanjanaporn Tansri, Sarawanee Choomklang, Pattrawut Chansangiam. Conjugate gradient algorithm for consistent generalized Sylvester-transpose matrix equations[J]. AIMS Mathematics, 2022, 7(4): 5386-5407. doi: 10.3934/math.2022299

    Related Papers:

    [1] Nunthakarn Boonruangkan, Pattrawut Chansangiam . Convergence analysis of a gradient iterative algorithm with optimal convergence factor for a generalized Sylvester-transpose matrix equation. AIMS Mathematics, 2021, 6(8): 8477-8496. doi: 10.3934/math.2021492
    [2] Fengxia Zhang, Ying Li, Jianli Zhao . The semi-tensor product method for special least squares solutions of the complex generalized Sylvester matrix equation. AIMS Mathematics, 2023, 8(3): 5200-5215. doi: 10.3934/math.2023261
    [3] Shousheng Zhu . Double iterative algorithm for solving different constrained solutions of multivariate quadratic matrix equations. AIMS Mathematics, 2022, 7(2): 1845-1855. doi: 10.3934/math.2022106
    [4] Huiling Wang, Nian-Ci Wu, Yufeng Nie . Two accelerated gradient-based iteration methods for solving the Sylvester matrix equation AX + XB = C. AIMS Mathematics, 2024, 9(12): 34734-34752. doi: 10.3934/math.20241654
    [5] Sani Aji, Poom Kumam, Aliyu Muhammed Awwal, Mahmoud Muhammad Yahaya, Kanokwan Sitthithakerngkiet . An efficient DY-type spectral conjugate gradient method for system of nonlinear monotone equations with application in signal recovery. AIMS Mathematics, 2021, 6(8): 8078-8106. doi: 10.3934/math.2021469
    [6] Habibu Abdullahi, A. K. Awasthi, Mohammed Yusuf Waziri, Issam A. R. Moghrabi, Abubakar Sani Halilu, Kabiru Ahmed, Sulaiman M. Ibrahim, Yau Balarabe Musa, Elissa M. Nadia . An improved convex constrained conjugate gradient descent method for nonlinear monotone equations with signal recovery applications. AIMS Mathematics, 2025, 10(4): 7941-7969. doi: 10.3934/math.2025365
    [7] Jamilu Sabi'u, Ibrahim Mohammed Sulaiman, P. Kaelo, Maulana Malik, Saadi Ahmad Kamaruddin . An optimal choice Dai-Liao conjugate gradient algorithm for unconstrained optimization and portfolio selection. AIMS Mathematics, 2024, 9(1): 642-664. doi: 10.3934/math.2024034
    [8] Sourav Shil, Hemant Kumar Nashine . Positive definite solution of non-linear matrix equations through fixed point technique. AIMS Mathematics, 2022, 7(4): 6259-6281. doi: 10.3934/math.2022348
    [9] Haixia Chang, Chunmei Li, Longsheng Liu . Generalized low-rank approximation to the symmetric positive semidefinite matrix. AIMS Mathematics, 2025, 10(4): 8022-8035. doi: 10.3934/math.2025368
    [10] Abdur Rehman, Ivan Kyrchei, Muhammad Zia Ur Rahman, Víctor Leiva, Cecilia Castro . Solvability and algorithm for Sylvester-type quaternion matrix equations with potential applications. AIMS Mathematics, 2024, 9(8): 19967-19996. doi: 10.3934/math.2024974
  • We develop an effective algorithm to find a well-approximate solution of a generalized Sylvester-transpose matrix equation where all coefficient matrices and an unknown matrix are rectangular. The algorithm aims to construct a finite sequence of approximated solutions from any given initial matrix. It turns out that the associated residual matrices are orthogonal, and thus, the desire solution comes out in the final step with a satisfactory error. We provide numerical experiments to show the capability and performance of the algorithm.



    Sylvester-type matrix equations show up naturally in several branches of mathematics and engineering. Indeed, many problems in vibration and structural analysis, robotics control and spacecraft control can be represented by the following general dynamical linear model:

    s1i=0Aix(i)+s2j=0Bju(j)=0 (1.1)

    where xRm×1 and uRn×1 are the state vector and the input vector respectively, and AiRm×m and BjRn×n are the system coefficient matrices; see e.g., [1,2]. The dynamical linear system (1.1) includes

    A1˙x+A0x+B0u=0,thedescriptorlinearsystem, (1.2)
    A2¨x+A1˙x+A0x+B0u=0,thesecondorderlinearsystem, (1.3)
    Akxk+Ak1xk1++A0x=Bu,thehighorderdynamicallinearsystem, (1.4)

    as special cases. Certain problems in control engineering, such as pole/eigenstructure assignment and observer design of the system (1.1) are closely related to the Lyapunov matrix equation AX+XAT=B, the Sylvester matrix equation AX+XB=C, and other realted equations; see e.g., [2,3,4,5,6,7]. In particular, the Sylvester matrix equation also plays a fundamental role in signal processing and model reduction; see e.g., [8,9,10]. These equations are special cases of a generalized Sylvester-transpose matrix equation:

    si=1AiXBi+tj=1CjXTDj=E. (1.5)

    A traditional way to solve Eq (1.5) for an exact solution is to transform it to an equivalent linear system via the Kronecker linearization; see Section 3 for details. However, this approach is only suitable when the dimensions of coefficient matrices are small. In practice, for a large-dimension case, it is enough to find a well-approximate solution via an iterative procedure, so that it is not necessary required memories as massive as traditional methods. There are several articles that consider problems that approximate the generalized Sylvester resistive matrix equations and constructs a finite sequence of approximated solutions from any given initial matrix. In the last five years, many researchers developed iterative methods for solving Sylvester-type matrix equations related to Eq (1.5). A group of Hermitian and skew-Hermitian splitting (HSS) methods aims to decompose a square matrix as the sum of its Hermitian part and skew-Hermitian part. There are several variantions of HSS, namely, GMHSS [11], preconditioned HSS [12], FPPSS [13], and ADSS [14]. A group of gradient-based iterative (GI) algorithms aims to construct a sequence of approximated solutions converging to the exact solution, based on the gradients of quadratic norm-error functions. The original GI method fora generalized Sylvester matrix equation was developed by Ding et al. [15]. Then Niu et al. [16] adjusted the GI method by introducing a weighted factor. After that a haif-step-update of GI method, called MGI method, introducing by Wang et al. [17]. The idea of GI algorithm can be used in conjuction with matrix diagonal-extraction to get AJGI [18] and MJGI [19] algorithms. See more GI algorithms for a generalized Sylvester matrix equations in [20,21,22]. For a generalized Sylvester-transpose matrix equation AXB+CXTD=E, there are GI algorithm [23] and an accelerated gradient-based iterative (AGBI) algorithm [24] to construct approximate solutions. There are also GI techniques based on optimization, e.g., [25,26,27]. See more computational methods for linear matrix equations in a survey [28]. The iterative procedures can be used to find solutions of certain nonlinear matrix equations; see e.g., [29,30,31]. There are also applications of such techniques to parameter estimation in dynamical systems; see e.g., [32,33,34].

    An idea of conjugate gradient (CG) is a remarkable technique constructing an orthogonal basis from the gradient of the associated quadratic function. There are several variations of CG to solve such matrix equations, e.g., BiCG [35], Bi-conjugate residual method [35], CGS [36], preconditioned nested splitting CG [37], generalized conjugate direction (GCD) method [38], conjugate gradient least-squares (CGLS) method [39], and GPBiCG [40].

    In this paper, we propose a conjugate gradient algorithm to solve the generalized Sylvester-transpose matrix Eq (1.5) in the consistent case, where all given coefficient matrices and the unknown matrix are rectangular (see Section 4). The algorithm aims to construct a sequence of approximate solutions of (1.5) from any given initial value. It turns out that a desire solution comes out in the final step of iterations with a satisfactory error (see Section 4). To validate the theory, we provide numerical experiments to show the applicability and the performance of the algorithm (see Section 5). In particular, the performance of the algorithm is significantly better than that of the direct Kronecker linearization and recent gradient-based iterative algorithms.

    In this section, we recall useful tools and facts from matrix analysis that are used in later discussions. Throughout, we denote the set of all m-by-n real matrices by Rm×n.

    Recall that the Kronecker product of A=[aij]Rm×n and BRp×q is defined by

    AB=[aijB]Rmq×np.

    Lemma 1 (see, e.g., [41]). The following properties hold for any compatible matrices A,B,C:

    1) (AB)T=ATBT,

    2) (A+B)C=AC+BC,

    3) A(B+C)=AB+AC.

    The vector operator Vec() assigns to each matrix A=[aij]Rm×n the column vector

    VecA=[a11am1a12am2a1namn]TRmn.

    This operator is bijective, linear, and compatible with the usual matrix multiplication in the following sense.

    Lemma 2 (see, e.g., [41]). For any ARm×n, BRp×q and XRn×p, we have

    VecAXB=(BTA)VecX.

    Recall that the commutation matrix P(m,n) is a permutation matrix defined by

    P(m,n)=mi=1nj=1EijETijRmn×mn (2.1)

    where each EijRm×n has entry 1 in (i,j)-th position and all other entries are 0.

    Lemma 3 (see, e.g., [41]). For any ARm×n and BRp×q, we have

    BA=P(n,p)T(AB)P(n,q). (2.2)

    Lemma 4 (see, e.g., [41]). For any matrix XRm×n, we have

    Vec(XT)=P(m,n)Vec(X). (2.3)

    Lemma 5 (see, e.g., [41]). For any matrices A,B,X,Y of compatible dimensions, we have

    (Vec(Y))T(AB)Vec(X)=tr(ATYTBX). (2.4)

    Recall that the Frobenius norm of ARm×n is defined by

    A=(mi=1nj=1a2ij)12=(tr(ATA))12.

    From now on, let m,n,p,q,r,s,k,N be such that mq=np. Consider the generalized Sylvester-transpose matrix Eq (1.5) where AiRm×n, BiRp×q, CjRm×p, DjRn×q, DRm×q, ERm×q are given matrices, and XRn×p is unknown. The Eq (1.5) includes the Lyapunov equation, the Sylvester equation, the equation AXB+CXD=E, and the equation AXB+CXTD=E as special cases.

    A direct method to solve Eq (1.5) is to transform it to an equivalent linear system. For convenience, denote P=P(n,p). By taking the vector operator to (1.5) and utilizing Lemma 4, we get

    VecE=Vec(si=1AiXBi+tj=1CjXTDj)=si=1(BTiAi)VecX+tj=1(DTjCj)VecXT=si=1(BTiAi)VecX+tj=1(DTjCj)PVecX=(si=1(BTiAi)+tj=1(DTjCj)P)VecX. (3.1)

    Let us denote x=VecX, b=VecE, and

    K=si=1(BTiAi)+tj=1(DTjCj)PRmq×np. (3.2)

    Thus, Eq (3.1) is equivalent to a linear algebraic system Kx=b. Hence, Eq (3.1) is consistent if and only if the associated linear system is consistent (i.e., rank[Kb]=rankK). When we solve for x, we can get the unknown matrix X using the injectivity of the vector operator. However, if the matrix coefficients Ai,Bi,Cj,Dj are of large sizes, then the size of K can be very large due to the Kronecker multiplication. Thus, traditional methods such as Gaussian elimination and LU factorization require a large memory to solve the linear system for an exact solution. Thus, the direct method is suitable for matrices of small sizes. For matrices of moderate/large sizes, it is enough to find a well-approximate solution for Eq (3.1) via an iterative procedure.

    The main task is to find a well-approximate solution of the matrix Eq (1.5):

    Problem 4.1. Let m,n,p,q,r,s,k,N be such that mq=np. Consider the generalized Sylvester-transpose matrix Eq (1.5) where AiRm×n, BiRp×q, CjRm×p, DjRn×q, DRm×q, ERm×q are given matrices, and XRn×p is unknown. Suppose that Eq (1.5) has a solution. Given an error ϵ>0, find ˜XRn×p such that

    Esi=1Ai˜XBitj=1Cj˜XTDj<ϵ.

    We will solve Problem 1 under an additional assumption that K in (3.2) is symmetric. We propose the following algorithm:

    Algorithm 1: A CG algorithm for a generalized Sylvester-transpose matrix equation
    AiRm×n, BiRp×q, CjRm×p, DjRn×q for any i=1,2,,s, j=1,2,,t and ERm×q;
    Given ϵ>0, set k=0, U0=0. Choose X0Rn×p
    R0=Esi=1AiX0Bitj=1CjX0Dj

     | Show Table
    DownLoad: CSV

    We call Xk the approximate solution at the k-th step. The main computation

    Xk+1=Xk+Rk2αk+1Uk+1,

    means that the next approximation Xk+1 is the sum between the current one Xk and the search direction Uk+1 with the step size Rk2/αk+1.

    Remark 4.2. The stopping rule of Algorithm 1 is based on the size of the residual matrix Rk. One can impose another stopping criterion besides Rkϵ, e.g., the norm of the difference Xk+1Xk between sucessive iterates is small enough.

    Remark 4.3. Let us discuss the complexity analysis for Algorithm 1. For convenience, suppose that all matrices in Eq (1.5) are of sizes n×n. Each step of the algorithm requires the matrix addition (O(n2)), the matrix multiplication (O(n3)), the matrix transposition (O(n2)), the matrix norm (O(n)), and the matrix trace (O(n)). In summary, the complexity analysis for each step is O(n3), and thus the algorithm runtime complexity is cubic time.

    Next, we will show that, for any given initial matrix X0, Algorithm 1 produces approximate solutions so that the set of residual matrices is an orthogonal set and, thus, we get the disire solution in a finite step. We divides the proof into several lemmas.

    Lemma 6. Consider Problem 1. Suppose that the sequence {Ri} is generated by Algorithm 1. We have

    Rk+1=RkRk2αk+1Vk+1fork=1,2,... (4.1)

    Proof. From Algorithm 1, we have that for any k,

    Rk+1=Esi=1AiXk+1Bitj=1CjXTk+1Dj=Esi=1Ai(Xk+Rk2αk+1Uk+1)Bitj=1Cj(XTk+Rk2αk+1UTk+1)Dj=Esi=1AiXkBitj=1CjXTkDjRk2αk+1(si=1AiUkBi+tj=1CjUTkDj)=RkRk2αk+1Vk+1.

    Lemma 7. Assume that the matrix K in (3.2) is symmetric. The sequences {Ui} and {Vi} generated by Algorithm 1 satisfy

    tr(UTmVn)=tr(VTmUn)foranym,n. (4.2)

    Proof. Applying Lemmas 1–5 and the symmetry of K, we have

    tr(VTmUn)=(VecVm)TVecUn=[si=1(BTiAi)VecUm+tj=1(DTjCj)VecUTm]TVecUn=[KVecUm]TVecUn=(VecUm)T KVecUn=(VecUm)T[Vec(si=1AiUnBi+tj=1CjUTnDj)]=(VecUm)TVecVn=tr(UTmVn)

    for any m,n.

    Lemma 8. Assume that the matrix K is symmetric. The sequences {Ri}, {Ui} and {Vi} are generated by Algorithm 1 satisfy

    tr(RTmRm1)=0andtr(UTm+1Vm)=0foranym. (4.3)

    Proof. To prove this conclusion, we use induction principle. In order to compute related terms, we use Lemmas 6 and 7. For m=1, we get

    tr(RT1R0)=tr[(R0R02α1V1)TR0]=tr(RT0R0)R02α1tr(VT1R0)=R02R02=0,

    and

    tr(UT2V1)=tr[(R1+R12R02U1)TV1]=tr(RT1V1)+R12R02tr(UT1V1)=α1R02tr(RT1R1)+α1R12R02=0.

    These imply that (4.3) hold for m=1.

    In the inductive step, for m=k we assume that tr(RTkRk1)=0 and tr(UTk+1Vk)=0. Then

    tr(RTk+1Rk)=tr[(RkRk2αk+1Vk+1)TRk]=tr(RTkRk)Rk2αk+1tr(VTk+1Rk)=tr(RTkRk)Rk2αk+1tr(VTk+1(Uk+1Rk2Rk12Uk))=Rk2Rk2αk+1tr(VTk+1Uk+1)=0,

    and

    tr(UTk+2Vk+1)=tr[(Rk+1+Rk+12Rk2Uk+1)TVk+1]=tr(RTk+1Vk+1)+Rk+12Rk2tr(UTk+1Vk+1)=tr(RTk+1(αk+1Rk2(Rk+1Rk)))+Rk+12Rk2αk+1=αk+1Rk2tr[(RTk+1Rk)(RTk+1Rk+1)]+Rk+12Rk2αk+1=0.

    Hence, Eq (4.3) holds for any m.

    Lemma 9. Assume that the matrix K is symmetric. Suppose the sequences {Ri}, {Ui} and {Vi} are generated by Algorithm 1.Then

    tr(RTmR0)=0,tr(UTm+1V1)=0foranym. (4.4)

    Proof. From Lemma 8 for m=1, we get tr(RT1R0)=0 and tr(PT2Q1)=0. Now, suppose that Eq (4.4) is true for all m=1,,k. From Lemmas 6 and 7, for m=k+1 we write

    tr(RTk+1R0)=tr[(RkRk2αk+1Vk+1)TR0]=tr(RTkR0)Rk2αk+1tr(VTk+1R0)=Rk2αk+1tr(VTk+1U1)=Rk2αk+1tr(UTk+1V1)=0,

    and

    tr(UTk+2V1)=tr(VTk+2U1)=tr[(αk+2Rk+12(Rk+2Rk+1))TU1]=αk+2Rk+12[tr(RTk+2U1)tr(RTk+1U1)]=αk+2Rk+12[tr(RTk+2R0)tr(RTk+1R0)]=0.

    Hence, Eq (4.4) holds for any m.

    Theorem 4.4. Assume that K is symmetric. Suppose the sequences {Ri}, {Ui} and {Vi} are generated by Algorithm 1. Then for any m,n such that mn, we have

    tr(RTm1Rn1)=0andtr(UTmVn)=0. (4.5)

    Proof. By Lemma 7 and the fact that tr(RTm1Rn1)=tr(RTn1Rm1) for any m,n, it suffices to prove (4.5) for any m,n such that m>n. By Lemma 8, Eq (4.5) holds for m=n+1. For m=n+2, we have

    tr(RTn+2Rn)=tr[(Rn+1Rn+12αn+2Vn+2)TRn]=Rn+12αn+2tr[VTn+2(Un+1Rn2Rn12Un)]=Rn+12αn+2Rn2Rn12[tr(Rn+1+Rn+12Rn2Un+1)TVn]=Rn+12αn+2Rn2Rn12[tr(αnRn12RTn+1(RnRn1))]=Rn+12αn+2Rn2Rn12αnRn12tr(RTn+1Rn1),
    tr(RTn+1Rn1)=tr[(RnRn2αn+1Vn+1)TRn1]=Rn2αn+1tr[VTn+1(UnRn12Rn22Un1)]=Rn2αn+1Rn12Rn22[tr(Rn+Rn2Rn12Un)TVn1]=Rn2αn+1Rn12Rn22[tr(αn1Rn22RTn(Rn1Rn2))]=Rn2αn+1Rn12Rn22αn1Rn22tr(RTnRn2),
    tr(UTn+2Vn)=tr[(Rn+1Rn+12Rn2Un+1)TVn]=tr[RTn+1(αnRn12(RnRn1))]=αnRn12tr[(RnRn2αn+1Vn+1)TRn1]=αnRn12Rn2αn+1[tr(VTn+1Un)Rn12Rn22tr(VTn+1Un1)]=αnRn12Rn2αn+1Rn12Rn22tr(UTn+1Vn1),

    and

    tr(UTn+1Vn1)=tr[(RnRn2Rn12Un)TVn1]=tr[RTn(αn1Rn22(Rn1Rn2))]=αn1Rn22tr[(Rn1Rn12αnVn)TRn2]=αn1Rn22Rn12αn[tr(VTnUn1)Rn22Rn32tr(VTnUn2)]=αn1Rn22Rn12αnRn22Rn32tr(UTnVn2),

    Similarly, we can write tr(RTn+2Rn) and tr(UTn+2Vn) in terms of tr(RTn+1Rn1) and tr(UTn+1Vn1), respectively. Repeating this process until the terms tr(RT2R0) and tr(UT3V1) show up. By Lemma 9, we get tr(RTn+2Rn)=0 and tr(UTn+2Vn)=0.

    Next, for m=n+3, we have

    tr(RTn+3Rn)=tr[(Rn+2Rn+22αn+3Vn+3)TRn]=Rn+22αn+3tr[VTn+3(Un+1Rn2Rn12Un)]=Rn+22αn+3Rn2Rn12[tr(Rn+2+Rn+22Rn+12Un+2)TVn]=Rn+22αn+3Rn2Rn12[tr(αnRn12RTn+2(RnRn1))]=Rn+22αn+3Rn2Rn12αnRn12tr(RTn+2Rn1),
    tr(RTn+2Rn1)=tr[(Rn+1Rn+12αn+2Vn+2)TRn1]=Rn+12αn+2tr[VTn+2(UnRn12Rn22Un1)]=Rn+12αn+2Rn12Rn22[tr(Rn+1+Rn+12Rn2Un+1)TVn1]=Rn+12αn+2Rn12Rn22[tr(αn1Rn22RTn+1(Rn1Rn2))]=Rn+12αn+2Rn12Rn22αn1Rn22tr(RTn+1Rn2),
    tr(UTn+3Vn)=tr[(Rn+2Rn+22Rn+12Un+2)TVn]=tr[RTn+2(αnRn12(RnRn1))]=αnRn12tr[(Rn+1Rn+12αn+2Vn+2)TRn1]=αnRn12Rn+12αn+2[tr(VTn+2Un)Rn12Rn22tr(VTn+2Un1)]=αnRn12Rn+12αn+2Rn12Rn22tr(UTn+2Vn1),

    and

    tr(UTn+2Vn1)=tr[(Rn+1Rn+12Rn2Un+1)TVn1]=tr[RTn+1(αn1Rn22(Rn1Rn2))]=αn1Rn22tr[(RnRn2αn+1Vn+1)TRn2]=αn1Rn22Rn2αn+1[tr(VTn+1Un1)Rn22Rn32tr(VTn+1Un2)]=αn1Rn22Rn2αn+1Rn22Rn32tr(UTn+1Vn2).

    Hence, we can write tr(RTn+3Rn) and tr(UTn+3Vn) in terms of tr(RTn+2Rn1) and tr(UTn+2Vn1), respectively. Repeating this process until the terms tr(RT3R0) and tr(U4V1) by Lemma 9, we get tr(RTn+3Rn)=0 and tr(UTn+3Vn)=0.

    Suppose that tr(RTm1RTn1)=tr(UTmVTn)=0 for m=n+1,,k. Then for m=k+1, we have

    tr(RTkRn1)=tr[(Rk1Rk12αkVk)TRn1]=tr(RTk1Rn1)Rk12αktr(VTkRn1)=Rk12αktr[VTk(UnRn12Rn22Un1)]=Rk12αk[tr(VTkUn)Rn12Rn22tr(VTkUn1)]=0.

    and

    tr(UTk+1Vn1)=tr[(Rk+Rk2Rk12Uk)TVn1]=tr(RTkVn1)+Rk2Rk12tr(UTkVn1)=tr[RTk(αn1Rn22(Rn1Rn2))]=αn1Rn22tr(RTkRn1RTkRn2)=0.

    Hence, tr(RTm1Rn1)=0 and tr(UTmVn)=0 for any m,n such that mn.

    Theorem 4.5. Consider Problem 4.1 under the assumption that the matrix K is symmetric. Suppose that the sequence {Xi} is generated by Algorithm 1. Then for given initial matrix X0Rn×p, an exact solution X can be obtained in at most np iteration steps.

    Proof. Suppose that Ri0 for i=0,1,,np1. Then we compute Xnp according to Algorithm 1. Assume that Rnp0. By Theorem 4.4, the set {R0,R1,...,Rnp} is orthogonal in Rn×p. So, {R0,R1,...,Rnp} is linearly independent. Since the dimension of Rn×p is np, any linearly independent subset of Rn×p must have at most np elements. So this is false because the set {R0,R1,...,Rnp} has np+1 elements. Thus, Rnp=0, hence Xnp is a solution of the equation.

    In this section, we report numerical results to illustrate the applicability and the effectiveness of Algorithm 1. All iterations have been carried out by MATLAB R2021a, on a macos (M1 chip 8C CPU/8C GPU/8GB/512GB). We perform the experiments for several generalized Sylvester-transpose matrix equations, and an interesting special case, namely, the Sylvester equation. We vary given coefficient matrices so that they are square/non-square sparse/dense matrices of moderate/large sizes. The dense matrices considered here are involved a matrix whose all entries are 1, which is denoted by ones. The identity matrix of size n×n is denoted by In. For each experiment, we set the stopping rule to be Rkϵ where ϵ=103. We discuss the performance of the algorithm through the norm of residual matrices, iteration number, and computational time (CT). The CT (in seconds) is measured by tic-toc function in MATLAB.

    In the following three examples, we concern the applicability of Algorithm 1 as well as its performance comparing to the direct Kronecker linearization mentioned in Section 3.

    Example 1. Consider a moderate-scaled generalized Sylvester-transpose equation

    A1XB1+A2XB2+C1XTD1+C2XTD2=E

    where all matrices are 50×50 tridiagonal matrices given by

    A1=tridiag(1,2,1),A2=tridiag(1,1,1),B1=tridiag(2,0,2),B2=tridiag(2,1,2),C1=tridiag(0,2,0),C2=tridiag(1,2,1),D1=tridiag(0,4,0),D2=tridiag(2,4,2),E=tridiag(1,1,9).

    We run Algorithm 1 using an initial matrix X0=0.25×onesR50×50. According to Theorem 4.5, Algorithm 1 will produce a solution of the equation within 104 iterations. The resulting simulation illustrated in Figure 1 shows the norms of residual matrices Rk at each iteration.

    Figure 1.  Relative error for Example 1.

    Althouh the errors Rk grow up and down during iterations, they generally climb down to zero. The algorithm takes 138 iterations to get a desire solution (so that Rk103), which is significantly less than the theoretical one (104 iterations). For the computational time, Algorithm 1 spends totally 0.131079 seconds to get a desire solution, while the direct Kronecker linearization consuming 1.581769 seconds to obtain the exact soluton. Thus, the performace of Algorithm 1 is significantly better than the direct method. Moreover, for sparse coefficient matrices, Agorithm 1 can produce a desire solution in a fewer iterations (that is, 138 iterations) than the theoretical one (that is, 104 iterations in this case).

    Example 2. Consider a generalized Sylvester-transpose matrix equation

    A1XB1+A2XB2+A3XB3+C1XTD1=E

    with rectangular coefficient matrices of moderate-scaled as follows:

    A1=tridiag(1,3,1),A2=tridiag(1,2,1),A3=tridiag(1,1,1)R40×40,B1=tridiag(2,1,2),B2=tridiag(1,3,1),B3=tridiag(2,3,2)R50×50,C1=3×onesR40×50,D1=3×onesR40×50,E=0.9×onesR40×50.

    Taking an initial matrix X0R40×50, we get an approximate solution XkR40×50 with a satisfactory error Rk103 in 164 steps, using 0.196250 seconds. We see in Figure 2 that during iterations, althouh the errors Rk grow up and down, they generally climb down to zero. On the other hand, the direct Kronecker linearization consumes 0.811170 seconds to get an exact solution. Thus, Agorithm 1 is applicable and effective.

    Figure 2.  Relative error for Example 2.

    Example 3. Consider a large-scaled generalized Sylvester-transpose equation

    A1XB1+C1XTD1+C2XTD2=E

    where all matrices are 100×100 tridiagonal matrices given by

    A1=tridiag(2,6,2),B1=tridiag(2,1,2),C1=tridiag(0,1,0),C2=tridiag(1,2,1),D1=tridiag(0,2,0),D2=tridiag(2,4,2),E=tridiag(1,8,1).

    The resulting simulation of Agorithm 1 using an initial matrix X0=0.5×onesR100×100 is shown in the next figure.

    Figure 3 shows the error gradually decreasing into ϵ=103 in 774 steps, consuming around 2 seconds.

    Figure 3.  Relative error for Example 3.

    Next, we investigate the effect of changing initial points. So we make experiments for the initial matrices X0=5×ones, X0=0, and X0=5×ones. Table 1 shows that, no matter the initial point, we get a desire solution in around 2 seconds. On the other hand the direct method consumes around 70 seconds to get an exacy solution. Thus, Agorithm 1 significantly outperforms the direct mathod.

    Table 1.  Relative error and CTs for Example 3.
    Initial matrix Iterations CT Relative error
    Direct - 69.500953 0
    X0=5×ones 830 2.247507 8.9145 ×104
    X0=0.5×ones 774 1.986782 7.5755 ×104
    X0=0 16 0.106482 5.3862 ×104
    X0=5×ones 830 2.190269 9.1232 ×104

     | Show Table
    DownLoad: CSV

    In the rest of numerical examples, we compare the performance of Algorithm 1 to the direct method as well as recent gradient-based iterative algorithms mentioned in Introduction.

    Example 4. Consider a large-scaled generalized Sylvester-transpose matrix equation

    AXB+CXTD=E,

    where A,B,C,D,E are 100×100 matrices as follows:

    A=tridiag(1,3,1),B=tridiag(1,7,1),C=6×ones,D=3×ones,E=0.7×I100.

    In fact, this equation has a unique solution. Despite the direct method, we compare the performance of Algorithm 1 to GI [23] and AGBI [24] algorithms. All iterative algorithms are implemented using the initial X0=0.001×I100R100×100.

    According to [23], the GI algorithm is applicable as long as a convergent factor μ satisfies

    0<μ<2λmax(AAT)λmax(BTB)+λmax(CCT)λmax(DTD),

    where λmax(AAT) is the largest eigenvalue of AAT. We run GI algorithm under 3 different convergent factors, namely, m1=6.1728×1012, m2=8.8183×1012 and m3=3.0864×1011. We implement AGBI algorithm with a convergent factor 0.000988 and a weighted factor 108. Figure 4 shows that the CG algorithm (Algorithm 1) converges faster than GI with m1, GI with m2, GI with m3, and AGBI algorithms. Table 2 shows that, in 30 iterations, GI, AGBI and the direct method consume a big amount of time to get the exact solution, while Algorithm 1 produces a small-error solution in a small time (0.073613 seconds).

    Figure 4.  Relative error for Example 4.
    Table 2.  Relative error and CTs for Example 4.
    Method Iterations CT Relative error
    CG 30 0.073613 0.000001
    GI with m1 30 0.109370 18.788266
    GI with m2 30 0.115890 16.853503
    GI with m3 30 0.109668 16.724640
    AGBI 30 0.118294 16.724926
    Direct - 66.928143 0

     | Show Table
    DownLoad: CSV

    Example 5. Consider a consistent generalized Sylvester-transpose matrix equation

    AXB+CXTD=E,

    with 100×100 coefficient matrices:

    A=tridiag(1,2,1),B=13×ones,C=3×ones,D=tridiag(3,6,3),E=1.2×ones.

    In fact, this matrix equation has a solution, which is not unique. We will seek for a solution of the equation using Algorithm 1, GI and AGBI algorithms with the same initial matrix X0=0.4×onesR100×100.

    We carry out GI algorithm with three different convergent factors, namely, m1=1.7132×108, m2=3.0837×108 and m3=1.5418×107. We implement AGBI algorithm with the convergent factor 0.000112 and the weighted factor 0.00005. Table 3 and Figure 5 express the computational time and the errors for 200 iterations of the simulations. We see that the computational time of CG algorithm is slightly less than those of GI (with parameters m1, m2, m3) and AGBI algorithms. However, the outcoming error produced by CG algorithm is significantly less than those of other algorithms.

    Table 3.  Relative error and CTs for Example 5.
    Method Iterations CT Relative error
    CG 200 0.395984 0.361597
    GI with m1 200 0.654457 1473.481117
    GI with m2 200 0.591405 1186.764341
    GI with m3 200 0.595052 645.799529
    AGBI 200 0.599059 1718.220885

     | Show Table
    DownLoad: CSV
    Figure 5.  Relative error for Example 5.

    Example 6. Consider the following Sylvester matrix equation

    AX+XB=C,

    where all coefficient matrices are 100×100 tridiagonal matrices given by

    A=tridiag(1,6,1),B=tridiag(3,0,3),C=tridiag(1,1,9).

    We compare the performance of CG algorithm (Algorithm 1) to GI [23], RGI [16], MGI [17] and AGBI [24] algorithms with parameters as shown in Table 4. We implement the algorithms with the same initial matrix X0=5×onesR100×100. Table 4 shows that the computational times for implementing 30 iterations of CG and other algorithms are close together. However, the relative errors in Figure 6 and Table 4 express that CG algorithm produces a sequence of well-approximate solutions in a few iterations with the lowest error comparing to other GI algorithms.

    Table 4.  Relative error and CTs for Example 6.
    Method Parameters Iterations CT Relative error
    Convergent factor weighted factor
    CG - - 10 0.018412 0.000000
    GI μ = 0.034482 - 10 0.019581 7.365534
    RGI μ = 0.266489 ω = 0.05 10 0.015338 8.786233
    MGI μ = 0.025316 - 10 0.019361 2.544848
    AGBI μ = 0.233918 ω = 0.05 10 0.022275 8.791699

     | Show Table
    DownLoad: CSV
    Figure 6.  Relative error for Example 6.

    We propose an iterative procedure (Algorithm 1) to construct a sequence of approximate solutions for the generalized Sylvester-transpose matrix Eq (1.5) with rectangular coefficient matrices. The algorithm is applicable whenever the matrix K, defined by Eq (3.2), is symmetric. In fact, the residual matrices Rk, produced by the algorithm, form an orthogonal set with respect to the usual inner product for matrices. Thus, we obtain the desire solution within a finite step, says, np steps. Numerical simulations have verified the applicability of the algorithm for square/non-square sparse/dense matrices of moderate/large sizes. The algorithm is always applicable no matter how we choose an initial matrix. Moreover, for sparse coefficient matrices of large size, the iteration number to get a desire solution can be dramatically less than np iterations. The performance of the algorithm is significantly better than the direct Kronecker linearization and recent gradient-based iterative algorithms when the matrix coefficients are of moderate/large sizes.

    This research project is supported by National Research Council of Thailand (NRCT): (N41A640234).

    All authors declare that they have no conflict of interest.



    [1] Y. Kim, H. S. Kim, J. Junkins, Eigenstructure assignment algorithm for second order systems, J. Guid. Control Dyn., 22 (1999), 729–731. http://dx.doi.org/10.2514/2.4444 doi: 10.2514/2.4444
    [2] B. Zhou, G. R. Duan, On the generalized Sylvester mapping and matrix equations, Syst. Control Lett., 57 (2008), 200–208. http://dx.doi.org/10.1016/j.sysconle.2007.08.010 doi: 10.1016/j.sysconle.2007.08.010
    [3] L. Dai, Singular control systems, Berlin: Springer, 1989.
    [4] G. R. Duan, Eigenstructure assignment in descriptor systems via output feedback: A new complete parametric approach, Int. J. Control., 72 (1999), 345–364. http://dx.doi.org/10.1080/002071799221154 doi: 10.1080/002071799221154
    [5] F. Lewis, A survey of linear singular systems, Circ. Syst. Signal Process., 5 (1986), 3–36. http://dx.doi.org/10.1007/BF01600184 doi: 10.1007/BF01600184
    [6] G. R. Duan, Parametric approaches for eigenstructure assignment in high-order linear systems, Int. J. Control Autom. Syst., 3 (2005), 419–429.
    [7] K. Nouri, S. Beik, L. Torkzadeh, D Baleanu, An iterative algorithm for robust simulation of the Sylvester matrix differential equations, Adv. Differ. Equ., 2020 (2020), http://dx.doi.org/10.1186/s13662-020-02757-z doi: 10.1186/s13662-020-02757-z
    [8] M. Epton, Methods for the solution of AXD - BXC = E and its applications in the numerical solution of implicit ordinary differential equations, BIT., 20 (1980), 341–345. http://dx.doi.org/10.1007/BF01932775 doi: 10.1007/BF01932775
    [9] D. Hyland, D. Bernstein, The optimal projection equations for fixed order dynamic compensation, IEEE Trans. Control., 29 (1984), 1034–1037. http://dx.doi.org/10.1109/TAC.1984.1103418 doi: 10.1109/TAC.1984.1103418
    [10] D. Calvetti, L. Reichel, Application of ADI iterative methods to the restoration of noisy images, SIAM J. Matrix Anal. Appl., 17 (1996), 165–186. http://dx.doi.org/10.1137/S0895479894273687 doi: 10.1137/S0895479894273687
    [11] M. Dehghan, A. Shirilord, A generalized modified Hermitian and skew-Hermitian splitting (GMHSS) method for solving complex Sylvester matrix equation, Appl. Math. Comput., 348 (2019), 632–651. http://dx.doi.org/10.1016/j.amc.2018.11.064 doi: 10.1016/j.amc.2018.11.064
    [12] S. Y. Li, H. L. Shen, X. H. Shao, PHSS Iterative method for solving generalized Lyapunov equations, Mathematics, 7 (2019), 38. http://dx.doi.org/10.3390/math7010038 doi: 10.3390/math7010038
    [13] H. L. Shen, Y. R. Li, X. H. Shao, The four-parameter PSS method for solving the Sylvester equation, Mathematics, 7 (2019), 105. http://dx.doi.org/10.3390/math7010105 doi: 10.3390/math7010105
    [14] M. Dehghan, A. Shirilord, Solving complex Sylvester matrix equation by accelerated double-step scale splitting (ADSS) method, Engineering with Computers, 37 (2021), 489–508. http://dx.doi.org/10.1007/s00366-019-00838-6 doi: 10.1007/s00366-019-00838-6
    [15] F. Ding, T. Chen, Gradient based iterative algorithms for solving a class of matrix equations, IEEE Trans. Automat. Comtr., 50 (2005), 1216–1221. http://dx.doi.org/10.1109/TAC.2005.852558 doi: 10.1109/TAC.2005.852558
    [16] Q. Niu, X. Wang, L. Z. Lu, A relaxed gradient based algorithm for solving Sylvester equation, Asian J. Control, 13 (2011), 461–464. http://dx.doi.org/10.1002/asjc.328 doi: 10.1002/asjc.328
    [17] X. Wang, L. Dai, D. Liao, A modified gradient based algorithm for solving Sylvester equation, Appl. Math. Comput., 218 (2012), 5620–5628. http://dx.doi.org/10.1016/j.amc.2011.11.055 doi: 10.1016/j.amc.2011.11.055
    [18] Z. Tian, M. Tian, C. Gu, X. Hao, An accelerated Jacobi-gradient based iterative algorithm for solving Sylvester matrix equations, Filomat, 31 (2017), 2381–2390. http://dx.doi.org/10.2298/FIL1708381T doi: 10.2298/FIL1708381T
    [19] N. Sasaki, P. Chansangiam, Modified Jacobi-gradient iterative method for generalized Sylvester matrix equation, Symmetry, 12 (2020), 1831. http://dx.doi.org/10.3390/sym12111831 doi: 10.3390/sym12111831
    [20] X. Zhang, X. Sheng, The relaxed gradient based iterative algorithm for the symmetric (skew symmetric) solution of the Sylvester equation AX + XB = C, Math. Probl. Eng., 2017 (2017), 1624969. http://dx.doi.org/10.1155/2017/1624969 doi: 10.1155/2017/1624969
    [21] A. Kittisopaporn, P. Chansangiam, W. Lewkeeratiyukul, Convergence analysis of gradient-based iterative algorithms for a class of rectangular Sylvester matrix equation based on Banach contraction principle, Adv. Differ. Equ., 2021 (2021), 17. http://dx.doi.org/10.1186/s13662-020-03185-9 doi: 10.1186/s13662-020-03185-9
    [22] N. Boonruangkan, P. Chansangiam, Convergence analysis of a gradient iterative algorithm with optimal convergence factor for a generalized Sylvester-transpose matrix equation, AIMS Mathematics, 6 (2021), 8477–8496. http://dx.doi.org/10.3934/math.2021492 doi: 10.3934/math.2021492
    [23] L. Xie, J. Ding, F. Ding, Gradient based iterative solutions for general linear matrix equations, Comput. Math. Appl., 58 (2009), 1441–1448. http://dx.doi.org/10.1016/j.camwa.2009.06.047 doi: 10.1016/j.camwa.2009.06.047
    [24] Y. J. Xie, C. F. Ma, The accelerated gradient based iterative algorithm for solving a class of generalized Sylvester ­transpose matrix equation, Appl. Math. Comp., 273 (2016), 1257–1269. http://dx.doi.org/10.1016/j.amc.2015.07.022 doi: 10.1016/j.amc.2015.07.022
    [25] A. Kittisopaporn, P. Chansangiam, Gradient-descent iterative algorithm for solving a class of linear matrix equations with applications to heat and Poisson equations, Adv. Differ. Equ., 2020 (2020), 324. http://dx.doi.org/10.1186/s13662-020-02785-9 doi: 10.1186/s13662-020-02785-9
    [26] A. Kittisopaporn, P. Chansangiam, The steepest descent of gradient-based iterative method for solving rectangular linear system with an application to Poisson's equation, Adv. Differ. Equ., 2020 (2020), 259. http://dx.doi.org/10.1186/s13662-020-02715-9 doi: 10.1186/s13662-020-02715-9
    [27] Y. Qi, L. Jin, H. Li, Y. Li, M. Liu, Discrete computational neural dynamics models for solving time-dependent Sylvester equation with applications to robotics and MIMO systems, IEEE Trans. Ind. Inform., 16 (2020), 6231–6241. http://dx.doi.org/10.1109/TII.2020.2966544 doi: 10.1109/TII.2020.2966544
    [28] V. Simoncini, Computational methods for linear matrix equations, SIAM Rev., 58 (2016), 377–441. http://dx.doi.org/10.1137/130912839 doi: 10.1137/130912839
    [29] H. Zhang, H. Yin, Refinements of the Hadamard and Cauchy Schwarz inequalities with two inequalities of the principal angles, J. Math. Inequal., 13 (2019), 423–435. http://dx.doi.org/10.7153/jmi-2019-13-28 doi: 10.7153/jmi-2019-13-28
    [30] H. Zhang, Quasi gradient-based inversion-free iterative algorithm for solving a class of the nonlinear matrix equations, Comput. Math. Appl., 77 (2019), 1233–1244. http://dx.doi.org/10.1016/j.camwa.2018.11.006 doi: 10.1016/j.camwa.2018.11.006
    [31] H. Zhang, L. Wan, Zeroing neural network methods for solving the Yang-Baxter-like matrix equation, Neurocomputing, 383 (2020), 409–418. http://dx.doi.org/10.1016/j.neucom.2019.11.101 doi: 10.1016/j.neucom.2019.11.101
    [32] F. Ding, G. Liu, X. Liu, Parameter estimation with scarce measurements, Automatica, 47 (2011), 1646–1655. http://dx.doi.org/10.1016/j.automatica.2011.05.007 doi: 10.1016/j.automatica.2011.05.007
    [33] F. Ding, Y. Liu, B. Bao, Gradient based and least squares based iterative estimation algorithms for multi-input multi-output systems, P. I. Mech. Eng. I-J. Sys., 226 (2012), 43–55. http://dx.doi.org/10.1177/0959651811409491 doi: 10.1177/0959651811409491
    [34] F. Ding, Combined state and least squares parameter estimation algorithms for dynamic systems, Appl. Math. Model., 38 (2014), 403–412. http://dx.doi.org/10.1016/j.apm.2013.06.007 doi: 10.1016/j.apm.2013.06.007
    [35] M. Hajarian, Developing Bi­CG and Bi­CR methods to solve generalized Sylvester-transpose matrix equations, Int. J. Autom. Comput., 11 (2014), 25–29. http://dx.doi.org/10.1007/s11633-014-0762-0 doi: 10.1007/s11633-014-0762-0
    [36] M. Hajarian, Matrix form of the CGS method for solving general coupled matrix equations, Appl. Math. Lett., 34 (2014), 37–42. http://dx.doi.org/10.1016/j.aml.2014.03.013 doi: 10.1016/j.aml.2014.03.013
    [37] Y. F. Ke, C. F. Ma, A preconditioned nested splitting conjugate gradient iterative method for the large sparse generalized Sylvester equation, Appl. Math. Comput., 68 (2014), 1409–1420. http://dx.doi.org/10.1016/j.camwa.2014.09.009 doi: 10.1016/j.camwa.2014.09.009
    [38] M. Hajarian, Generalized conjugate direction algorithm for solving the general coupled matrix equations over symmetric matrices, Numer. Algor., 73 (2016), 591–609. http://dx.doi.org/10.1007/s11075-016-0109-8 doi: 10.1007/s11075-016-0109-8
    [39] M. Hajarian, Extending the CGLS algorithm for least squares solutions of the generalized Sylvester-transpose matrix equations, J. Franklin Inst., 353 (2016), 1168–1185. http://dx.doi.org/10.1016/j.jfranklin.2015.05.024 doi: 10.1016/j.jfranklin.2015.05.024
    [40] M. Dehghan, R. Mohammadi-Arani, Generalized product-type methods based on Bi-conjugate gradient (GPBiCG) for solving shifted linear systems, Comput. Appl. Math., 36 (2017), 1591–1606. http://dx.doi.org/10.1007/s40314-016-0315-y doi: 10.1007/s40314-016-0315-y
    [41] R. Horn, C. Johnson, Topics in matrix analysis, Cambridge: Cambridge University Press, 1991. http://dx.doi.org/10.1017/CBO9780511840371
  • This article has been cited by:

    1. Kanjanaporn Tansri, Pattrawut Chansangiam, Conjugate Gradient Algorithm for Least-Squares Solutions of a Generalized Sylvester-Transpose Matrix Equation, 2022, 14, 2073-8994, 1868, 10.3390/sym14091868
    2. Victor G. Lopez, Matthias A. Müller, 2023, An Efficient Off-Policy Reinforcement Learning Algorithm for the Continuous-Time LQR Problem, 979-8-3503-0124-3, 13, 10.1109/CDC49753.2023.10384256
    3. Rui Qi, Caiqin Song, Iterative algorithm to the generalized coupled Sylvester‐transpose matrix equations with application in robust and minimum norm observer design of linear systems, 2024, 45, 0143-2087, 2008, 10.1002/oca.3134
    4. Janthip Jaiprasert, Pattrawut Chansangiam, Exact and least-squares solutions of a generalized Sylvester-transpose matrix equation over generalized quaternions, 2024, 32, 2688-1594, 2789, 10.3934/era.2024126
  • Reader Comments
  • © 2022 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(2709) PDF downloads(107) Cited by(4)

Figures and Tables

Figures(6)  /  Tables(4)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog