# LSMR Iterative Method for General Coupled Matrix Equations

Authors: F. Toutounian; D. Khojasteh Salkuyeh; M. Mojarrab
Journal: Journal of Applied Mathematics
 (2015)

---

## Abstract

By extending the idea of LSMR method, we present an iterative method to solve the general coupled matrix equations∑k=1qAikXkBik=Ci, i=1,2,…,p, (including the generalized (coupled) Lyapunov and Sylvester matrix equations as special cases) over some constrained matrix groups (X1,X2,…,Xq), such as symmetric, generalized bisymmetric, and (R,S)-symmetric matrix groups. By this iterative method, for any initial matrix group  (X1(0),X2(0),…,Xq(0)), a solution group (X1*,X2*,…,Xq*) can be obtained within finite iteration steps in absence of round-off errors, and the minimum Frobenius norm solution or the minimum Frobenius norm least-squares solution group can be derived when an appropriate initial iterative matrix group is chosen. In addition, the optimal approximation solution group to a given matrix group (X¯1,X¯2,…,X¯q) in the Frobenius norm can be obtained by finding the least Frobenius norm solution group of new general coupled matrix equations. Finally, numerical examples are given to illustrate the effectiveness of the presented method.


---

## Body

## 1. Introduction

In control and system theory [1–7], we often encounter Lyapunov and Sylvester matrix equations which have been playing a fundamental role. Due to the important roles of the matrix equations, the studies on the matrix equations have been addressed in a large body of papers [8–27]. By using the hierarchical identification principle [9–11, 28–32], a gradient-based iterative (GI) method was presented to solve the solutions and the least-squares solutions of the general coupled matrix equations. In [19, 33], Zhou et al. deduced the optimal parameter of the GI method for solving the solutions and the weighted least-squares solutions of the general coupled matrix equations. Dehghan and Hajarian [34–36] introduced several iterative methods to solve various linear matrix equations.In [12, 17], Huang et al. presented finite iterative algorithms for solving generalized coupled Sylvester systems. Li and Huang [37] proposed a matrix LSQR iterative method to solve the constrained solutions of the generalized coupled Sylvester matrix equations. Hajarian [38] presented the generalized QMRCGSTAB algorithm for solving Sylvester-transpose matrix equations. Recently, Lin and Simoncini [39] established minimal residual methods for large scale Lyapunov equations. They explored the numerical solution of this class of linear matrix equations when a minimal residual (MR) condition is used during the projection step.In this paper, we construct a matrix iterative method based on the LSMR algorithm [40] to solve the constrained solutions of the following problemsCompatible matrix equations are as follows:(1)∑k=1q‍AikXkBik=Ci,i=1,2,…,p.Least-squares problem is as follows:(2)min⁡∑i=1p‍∑k=1q‍AikXkBik-Ci21/2.Matrix nearness problem is as follows:(3)min⁡X1,X2,…,Xq∈SE∑k=1q‍Xk-X¯k2,where Aik, Bik, Ci, i=1,2,…,p, k=1,2,…,q, are constant matrices with suitable dimensions, Xk, k=1,2,…,q, are unknown matrices to be solved, X¯k, k=1,2,…,q, are given matrices, and SE is the solution set of (1) or problem (2).This paper is organized as follows. In Section2, we will briefly review the LSMR algorithm for solving linear systems of equations. In Section 3, we propose the matrix LSMR iterative algorithms for solving the problems (1)-(2). In Section 4, we solve the problem (3) by finding the minimum Frobenius norm solution group of the corresponding new general coupled matrix equations. In Section 5, numerical examples are given to illustrate the efficiency of the proposed iterative method. Finally, we make some concluding remarks in Section 6.The notations used in this paper can be summarized as follows.tr⁡(A) represents the trace of the matrix A. For A,B∈Rm×n, notation A⊗B is Kronecker product and 〈A,B〉=tr⁡(BTA) is the inner product with the Frobenius norm A=〈A,A〉=tr⁡(ATA). The use of vec⁡(A) represents the vector operator defined as(4)vec⁡(A)=a1Ta2T⋯anTT,where ak is the kth column of A. The generalized bisymmetric matrices, the (R,S)-symmetric matrices, and the symmetric orthogonal matrices can be defined as follows.Definition 1 (see [41]).
A matrixP∈Rn×n is said to be a symmetric orthogonal matrix (P∈SORn×n), if PT=P and P2=In.Definition 2 (see [42]).
For given symmetric orthogonal matricesR∈Rm×m, S∈Rn×n, we say a matrix X∈Rm×n is (R,S)-symmetric (X∈RSSRn×n), if X=RXS.Definition 3 (see [43]).
For a given symmetric orthogonal matrixP∈Rn×n, a matrix X∈Rn×n is said to be a generalized bisymmetric matrix (X∈GBRSn×n), if XT=X and X=PXP.


## 2. LSMR Algorithm

In this section, we briefly review some fundamental properties of the LSMR algorithm [40], which is an iterative method for computing a solution x to either of the following problems.Compatible linear systems are as follows:(5)Mx=f.Least-squares problem is as follows:(6)min⁡Mx-f2,where M∈Rm×n and f∈Rm. The LSMR algorithm uses an algorithm of Golub and Kahan [44], which stated as procedure Bidiag 4, to reduce M to lower bidiagonal form. The procedure Bidiag 4 can be described as follows.Bidiag 4 (starting vectorf; reduction to lower bidiagonal form).
(7)
β
1
u
1
=
f
,
α
1
v
1
=
M
T
u
1
,
β
i
+
1
u
i
+
1
=
M
v
i
-
α
i
u
i
,
α
i
+
1
v
i
+
1
=
M
T
u
i
+
1
-
β
i
+
1
v
i
,
m
m
m
m
m
m
m
m
i
=
1,2
,
…
.The scalers αi≥0 and βi≥0 are chosen such that ui2=vi2=1.
The following properties presented in [7] illustrate that the procedure Bidiag 4 has finite termination property.Property 1.
Suppose thatk steps of the procedure Bidiag 4 have been taken; then the vectors v1,v2,…,vk and u1,u2,…,uk,uk+1 are the orthonormal basis of the Krylov subspaces Kk(MTM,v1) and Kk+1(MMT,u1), respectively.Property 2.
The procedure Bidiag4 will stop at step m if and only if min⁡{μ,λ} is m, where μ is the grade of v1 with respect to MTM and λ is the grade of u1 with respect to MMT.By using the procedure Bidiag4, the LSMR method constructs an approximation solution of the form xk=Vkyk, where Vk=(v1,v2,…,vk), which solves the least-squares problem, min⁡ykMTrk, where rk=f-Mxk is the residual for the approximate solution xk. The main steps of the LSMR algorithm can be summarized as shown in Algorithm 1.Algorithm 1:LSMR algorithm.
Setβ1u1=f, α1v1=MTu1, α¯1=α1, ζ1¯=α1β1, ρ0=1, ρ¯0=1, c0¯=1, s0¯=0, h1=v1, h¯0=0, x0=0
Fork=1,2,…, until convergence Do:
βk+1uk+1=Mvk-αkuk
αk+1vk+1=MTuk+1-βk+1vk
ρk=(α¯k2+βk+12)1/2
ck=α¯k/ρk
sk=βk+1/ρk
θk+1=skαk+1
α¯k+1=ckαk+1
θ¯k=s¯k-1ρk
ρ¯k=((c¯k-1ρk)2+θk+12)1/2
c¯k=c¯k-1ρk/ρ¯k
s¯k=θk+1/ρ¯k
ζk=c¯kζ¯k
ζ¯k+1=-s¯kζ¯k
h¯k=hk-(θ¯kρk/(ρk-1ρ¯k-1))h¯k-1
xk=xk-1+(ζk/(ρkρ¯k))h¯k
hk+1=vk+1-(θk+1/ρk)hk
If |ζ¯k+1| is small enough then stop
End Do.More details about the LSMR algorithm can be found in [40].The stopping criterion may be used asMTrk=MTf-Mxk2=ζ¯k+1 for the compatible linear systems (5) and for the least-squares problem (6). Other stopping criteria can also be used and are not listed here. Reader can see [40] for details. Clearly, the sequence xk∈Range(MT) generated by the LSMR algorithm converges to the unique minimum norm solution of (5) or the unique minimum norm least-squares solution of problem (6).


## 3. A Matrix LSMR Iterative Method

In this section, we will present our matrix iterative method based on the LSMR algorithm, for solving (1) and problem (2). For the unknown matrices Xi∈Rn×n, i=1,2,…,q, by using the Kronecker product, (1) and problem (2) are equivalent to (5) and problem (6), respectively, with(8)M=B11T⊗A11B12T⊗A12⋯B1qT⊗A1qB21T⊗A21B22T⊗A22⋯B2qT⊗A2q⋮⋮⋮Bp1T⊗Ap1Bp2T⊗Ap2⋯BpqT⊗Apq,x=vec⁡(X1)vec⁡(X2)⋮vec⁡(Xq),f=vec⁡(C1)vec⁡(C2)⋮vec⁡(Cp).Hence, by using the invariance of the Frobenius norm under unitary transformations, it is easy to prove that the vector formβ1u1=f, α1v1=MTu1, βi+1ui+1=Mvi-αiui, and αi+1vi+1=MTui+1-βi+1vi (i=1,2,…) in LSMR algorithm can be rewritten in matrix forms, respectively, as(9)β1U1j=Cj,j=1,2,…,p,β1=∑j=1p‍Cj21/2,α1V1j=∑i=1p‍AijTU1iBijT,j=1,2,…,q,α1=∑k=1q‍∑i=1p‍AikTU1iBikT21/2,βi+1Ui+1j=∑k=1q‍AjkVikBjk-αiUij,j=1,2,…,p,βi+1=∑j=1p‍∑k=1q‍AjkVikBjk-αiUij21/2,αi+1Vi+1k=∑j=1p‍AjkTUi+1jBjkT-βi+1Vik,k=1,2,…,q,αi+1=∑k=1q‍∑j=1p‍AjkTUi+1jBjkT-βi+1Vik21/2.If the unknown matrices X1,X2,…,Xq∈GBRSn×n, then (1) and problem (2) are equivalent to (5) and problem (6), respectively, with(10)M=B11T⊗A11B12T⊗A12⋯B1qT⊗A1qB21T⊗A11B22T⊗A22⋯B2qT⊗A2q⋮⋮⋮Bp1T⊗Ap1Bp2T⊗Ap2⋯BpqT⊗ApqA11⊗B11TA12⊗B12T⋯A1q⊗B1qTA21⊗B21TA22⊗B22T⋯A2q⊗B2qT⋮⋮⋮Ap1⊗Bp1TAp2⊗Bp2T⋯Apq⊗BpqT(B11TP1)⊗(A11P1)(B12TP2)⊗(A12P2)⋯(B1qTPq)⊗(A1qPq)(B21TP1)⊗(A21P1)(B22TP2)⊗(A22P2)⋯(B2qTPq)⊗(A2qPq)⋮⋮⋮(Bp1TP1)⊗(Ap1P1)(Bp2TP2)⊗(Ap2P2)⋯(BpqTPq)⊗(ApqPq)(A11P1)⊗(B11TP1)(A12P2)⊗(B12TP2)⋯(A1qPq)⊗(B1qTPq)(A21P1)⊗(B21TP1)(A22P2)⊗(B22TP2)⋯(A2qPq)⊗(B2qTPq)⋮⋮⋮(Ap1P1)⊗(Bp1TP1)(Ap2P2)⊗(Bp2TP2)⋯(ApqPq)⊗(BpqTPq),x=vec⁡(X1)vec⁡(X2)⋮vec⁡(Xq),f=vec⁡(C1)vec⁡(C2)⋮vec⁡(Cp)vec⁡(C1T)vec⁡(C2T)⋮vec⁡(CpT)vec⁡(C1)vec⁡(C2)⋮vec⁡(Cp)vec⁡(C1T)vec⁡(C2T)⋮vec⁡(CpT),where P1,P2,…,Pq∈SORn×n. Hence, the vector forms of β1u1=f, α1v1=MTu1, βi+1ui+1=Mvi-αiui, and αi+1vi+1=MTui+1-βi+1vi (i=1,2,…) in LSMR algorithm can be rewritten in matrix forms, respectively, as(11)β1U1j=Cj,j=1,2,…,p,β1=2∑j=1p‍Cj21/2,α1V1(j)=∑i=1p‍AijTU1(i)BijT+BijU1(i)TAij+(PjAijT)U1(i)(BijTPj)+(PjBij)U1(i)T(AijPi),j=1,2,…,q,α1=∑k=1q‍∑i=1p‍AikTU1(i)BikT+BikU1(i)TAik+(PkAikT)U1(i)(BikTPk)+PkBikU1iTAikPi∑i=1p‍2∑i=1p‍1/2,βi+1Ui+1j=∑k=1q‍AjkVikBjk-αiUij,j=1,2,…,p,βi+1=2∑j=1p‍∑k=1q‍AjkVikBjk-αiUij21/2,αi+1Vi+1k=∑j=1p‍AjkTUi+1jBjkT+BjkUi+1jTAjk+PkAjkT×Ui+1jBjkTPk+PkBjkUi+1jTAjkPk-βi+1Vik,k=1,2,…,q,αi+1=∑k=1q‍∑j=1p‍AjkTUi+1jBjkT+BjkUi+1jTAjk+PkAjkTUi+1jBjkTPk+PkBjkUi+1jTAjkPk-βi+1Vik∑k=1q‍2∑k=1q‍1/2.If the unknown matrices X1,X2,…,Xq∈RSSRn×n, then (1) and problem (2) are equivalent to (5) and problem (6), respectively, with(12)M=B11T⊗A11B12T⊗A12⋯B1qT⊗A1qB21T⊗A21B22T⊗A22⋯B2qT⊗A2q⋮⋮⋮Bp1T⊗Ap1Bp2T⊗Ap2⋯BpqT⊗Apq(B11TS1)⊗(A11R1)(B12TS2)⊗(A12R2)⋯(B1qTSq)⊗(A1qRq)(B21TS1)⊗(A21R1)(B22TS2)⊗(A22R2)⋯(B2qTSq)⊗(A2qRq)⋮⋮⋮(Bp1TS1)⊗(Ap1R1)(Bp2TS2)⊗(Ap2R2)⋯(BpqTSq)⊗(ApqRq),x=vec⁡(X1)vec⁡(X2)⋮vec⁡(Xq),f=vec⁡(C1)vec⁡(C2)⋮vec⁡(Cp)vec⁡(C1)vec⁡(C2)⋮vec⁡(Cp),where Pi, Ri∈SORn×n,i=1,2,…,q. Hence, the vector forms of β1u1=f, α1v1=MTu1,  βi+1ui+1=Mvi-αiui, and αi+1vi+1=MTui+1-βi+1vi (i=1,2,…) in LSMR algorithm can be rewritten in matrix forms, respectively, as(13)β1U1j=Cj,j=1,2,…,p,β1=2∑j=1p‍Cj21/2,α1V1j=∑i=1p‍AijTU1iBijT+RjAijTU1iBijTSj,nj=1,2,…,q,α1=∑k=1q‍∑i=1p‍AikTU1iBikT+RkAikTU1iBikTSk21/2,βi+1Ui+1j=∑k=1q‍AjkVikBjk-αiUij,j=1,2,…,p,βi+1=2∑j=1p‍∑k=1q‍AjkVikBjk-αiUij21/2,αi+1Vi+1k=∑j=1p‍AjkTUi+1jBjkT+RkAjkTUi+1jBjkTSk-βi+1Vik,k=1,2,…,q,αi+1=∑k=1q‍∑j=1p‍(AjkTUi+1(j)BjkT+Rk(AjkTUi+1(j)BjkT)Sk)-βi+1Vik∑k=1q‍2∑k=1q‍1/2.If the unknown matricesX1,X2,…,Xq∈SRn×n (the set of symmetric matrices), then (1) and problem (2) are equivalent to (5) and problem (6), respectively, with(14)M=B11T⊗A11B12T⊗A12⋯B1qT⊗A1qB21T⊗A21B22T⊗A22⋯B2qT⊗A2q⋮⋮⋮Bp1T⊗Ap1Bp2T⊗Ap2⋯BpqT⊗ApqA11⊗B11TA12⊗B12T⋯A1q⊗B1qTA21⊗B21TA22⊗B22T⋯A2q⊗B2qT⋮⋮⋮Ap1⊗Bp1TAp2⊗Bp2T⋯Apq⊗BpqT,(15)x=vec⁡(X1)vec⁡(X2)⋮vec⁡(Xq),f=vec⁡(C1)vec⁡(C2)⋮vec⁡(Cp)vec⁡(C1T)vec⁡(C2T)⋮vec⁡(CpT).Hence, the vector forms of β1u1=f, α1v1=MTu1, βi+1ui+1=Mvi-αiui, and αi+1vi+1=MTui+1-βi+1vi (i=1,2,…) in LSMR algorithm can be rewritten in matrix forms, respectively, as(16)β1U1j=Cj,j=1,2,…,p,β1=2∑j=1p‍Cj21/2,α1V1j=∑i=1p‍AijTU1iBijT+BijU1iTAij,j=1,2,…,q,α1=∑k=1q‍∑i=1p‍AikTU1iBikT+BikU1iTAik21/2,βi+1Ui+1j=∑k=1q‍AjkVikBjk-αiUij,j=1,2,…,p,βi+1=2∑j=1p‍∑k=1q‍AjkVikBjk-αiUij21/2,αi+1Vi+1k=∑j=1p‍AjkTUi+1jBjkT+BjkUi+1jTAjk-βi+1Vik,mmmmmmmmmmmmmmk=1,2,…,q,αi+1=∑k=1q‍∑j=1p‍AjkTUi+1jBjkT+BjkUi+1jTAjk-βi+1Vik21/2.From above results, we can obtain the matrix form iteration method of LSMR algorithm for solving the constrained solution group of (1) and problem (2). When the unknown matrices X1,X2,…,Xq∈SRn×n, the matrix form iterative method is given as shown in Algorithm 2.Algorithm 2:LSMR_SR_M algorithm.
SetX0(k)=0, k=1,2,…,q
β1=2(∑j=1p‍‖Cj‖2)1/2, U1(j)=Cj/β1, j=1,2,…,p
α1=(∑k=1q‍‖∑i=1p‍AikTU1(i)BikT+BikU1(i)TAik‖2)1/2
V1(k)=∑i=1p‍(AikTU1(i)BikT+BikU1(i)TAik)/α1, k=1,2,…,q
Setα¯1=α1, ζ¯1=α1β1, ρ0=1, ρ¯0=1, c¯0=1, s¯0=0, h1(k)=V1(k),k=1,2,…,q, h¯0(k)=0,k=1,2,…,q
Fork=1,2,…, until converges Do:
βk+1=2(∑i=1p‍‖∑j=1q‍AijVk(j)Bij-αkUk(i)‖2)1/2
Uk+1(i)=∑j=1q‍AijVk(j)Bij-αkUk(i)/βk+1, i=1,2,…,p
αk+1=(∑j=1q‍‖∑i=1p‍(AijTUk+1(i)BijT+BijUk+1(i)TAij)-βk+1Vk(j)‖2)1/2
Vk+1(i)=(∑i=1p‍(AijTUk+1(i)BijT+BijUk+1(i)TAij)-βk+1Vk(j))/αk+1, i=1,2,…,q
ρk=(α¯k2+βk+12)1/2
ck=α¯k/ρk
sk=βk+1/ρk
θk+1=skαk+1
α¯k+1=ckαk+1
θ¯k=s¯k-1ρk
ρ¯k=((c¯k-1ρk)2+θk+12)1/2
c¯k=c¯k-1ρk/ρ¯k
s¯k=θk+1/ρ¯k
ζk=c¯kζ¯k
ζ¯k+1=-s¯kζ¯k
h¯k(i)=hk(i)-(θ¯kρk/(ρk-1ρ¯k-1))h¯k-1(i), i=1,2,…,q
Xk(i)=Xk-1(i)+(ζk/(ρkρ¯k))h¯k(i), i=1,2,…,q
hk+1(i)=Vk+1(i)-(θk+1/ρk)hk(i), i=1,2,…,q
If |ζ¯k+1| is small enough then stop
End Do.


## 4. The Solution Group of Problem (3)

Now, we consider the solution group of the matrix nearness problem (3) for given matrix group (X¯1,X¯2,…,X¯q), where X¯k∈Rn×n,  k=1,2,…,q. If Xk∈SRn×n, it is easy to prove that(17)min⁡Xk∈SRn×n,k=1,2,…,q∑k=1q‍Xk-X¯k2=min⁡Xk∈SRn×n,k=1,2,…,q∑k=1q‍Xk-X¯k+X¯kT22+∑k=1q‍X¯k-X¯kT22.Let(18)X^k=Xk-X¯k+X¯kT2,k=1,2,…,q,C^j=Cj-∑k=1q‍AjkX¯k+X¯kT2Bjk,j=1,2,…,p;then problem (3) is equivalent to finding the minimum Frobenius norm symmetric solution group or minimum Frobenius norm least-squares symmetric solution group of the following problems, respectively.Compatible matrix equations are as follows:(19)∑k=1q‍AikX^kBik=C^i,i=1,2,…,p.Least-squares problem is as follows:(20)min⁡∑i=1p‍∑k=1q‍AikX^kBik-C^i21/2.By LSMR_SR_M method, we can get the minimum Frobenius norm symmetric solution group (X^1*,X^2*,…,X^q*) of (19) (or the minimum Frobenius norm least-squares symmetric solution group of problem (20)). Then, the optimal approximate solution group (X~1,X~2,…,X~q) of problem (3) can be obtained; that is, X~k=X^k*+(X¯k+X¯kT)/2.


## 5. Numerical Examples

To compare the behavior of the proposed matrix method discussed in the previous section with the CGNE method [43] and the matrix LSQR iterative method (LSQR_M) [37], we present in this section numerical results for three examples. All the numerical computations are performed in MATLAB 7.Example 1.
Suppose that the matricesAij,  Bij,  Cj,  i,j=1,2, are given by the following matrices:(21)A11=8-12-345022103-724-2-1-12441121,B11=9-2-311054-2-225-4-8190251-3-3119,A12=4-34-41-2-26-449605154533917-59,B12=-2-1-2-32154969-3111-525414-8221,A21=123-52461-3944-2232281-9314-38,B21=923-1-13-1885-1-2-8-5-2-1-12731-1-882,A22=-92321-33943826-8812-94-22-1-2-2-3,B22=-1-29891-3-3-39-1-3-9-591-9-65-12-7-458.The C1 and C2 matrices are chosen such that X1=In and X2=En, where In and En are the n×n identity matrix and the n×n matrix whose entries are all one, respectively.In Figure1, we display the convergence curves of the function log⁡10δk, with(22)δk=max⁡R1kFR10F,R2kFR20F,where Ri(k),  i,k=1,2, is the residual matrix of the ith equation in kth iteration. The initial iterative matrices in all the iterative methods are chosen as zero matrices of suitable size. Figure 1 confirms that the proposed algorithm has faster convergence rate and higher accuracy than the CGNE method and similar behavior to the matrix LSQR iterative method.Figure 1
Convergence history of the LSMR, CGNE, and LSQR iterative methods for Example1.Example 2.
Suppose that the matricesAij,  Bij,  Cj,  i,j=1,2, are given by the following matrices:(23)A11=tridiag(-1,6,-1),B11=tridiag(1,8,-1),A12=(0.1)In,B12=tridiag(1,0,1),A21=0.1In,B21=tridiag-2,1,-2,A22=tridiag(-1,-3,-1),B22=tridiag(1,6,2).As Example 1, the C1 and C2 matrices are chosen such that X1=In and X2=En with n=400 and the initial iterative matrices in all the iterative methods are chosen as zero matrices of suitable size. In Figure 2, as Figure 1, we display the convergence curves of the function log⁡10δk,. This figure shows that the LSMR method outperforms the CGNE and LSQR methods.Figure 2
Convergence history of the LSMR, CGNE, and LSQR iterative methods for Example2.Example 3 (see [45]).
Consider the convection diffusion equation with the Dirichlet boundary conditions(24)Lu:=-Δu+2ν∂u∂x+2ν∂u∂y=fon⁡Ω,u=gon⁡∂Ω.Here Ω is the unit square [0,1]×[0,1]. The operator L was discretized using central finite differences on Ω, with mesh size h=1/(n+1) in the “x” direction and k=1/(p+1) in the “y” direction. This yields a linear system of algebraic equations that can be written as a Sylvester matrix equation(25)AX-XB=C,(as a particular case of (1) with A11=A,  A12=-In,  B11=Ip,  B12=B,  C1=C, and X1=X2=X) where tridiagonal matrices A and B are given by(26)A=-1h2tridiag1+νh,-2,1-νh,B=1k2tridiag1+νk,-2,1-νk.The right-hand side matrix C is obtained as follows:(27)Ci,j=fxi+1,yj+1,fori=2,…,n-1,mmmmmmmmmmj=2,…,p-1,C1,1=fx2,y2+1+νhh2g0,y2+1+νkk2gx2,0,C1,p=fx2,yp+1+1+νhh2g0,yp+1+1-νkk2gx2,1,Cn,1=fxn+1,y2+1-νhh2g1,y2+1+νkk2gxn+1,0,Cn,p=fxn+1,yp-1+1-νhh2g1,yp+1+1-νkk2gxn+1,1,C1,j=fx2,yj+1+1+νhh2g0,yj+1,mmmmmmforj=2,…,p-1,Cn,j=fxn+1,yj+1+1-νhh2g1,yj+1,mmmmmmmforj=2,…,p-1,Ci,1=fxi+1,y2+1+νkk2gxi+1,0,mmmmmmfori=2,…,n-1,Ci,p=f(xi+1,yp+1)+1-νkk2g(xi+1,1),mmmmmmmfori=2,…,n-1.In this example, the functions f and g were chosen such that the exact solution is (28)u(x,y)=xe-x2-y2on the domain Ω. In addition, we used the symmetric successive overrelaxation (SSOR) preconditioner for the matrix equation (25) to increase the convergence rate. It is easy to prove that the matrix equation (25) is equivalent to the np×np linear system:(29)A~x~=c~,where A~=Ip⊗A-BT⊗In,  c~=vec⁡(C), and x~=vec⁡(X).
The matricesA and B can be written as(30)A=DA-EA-FA,B=DB-EB-FB,where DA is the diagonal of A and -EA and -FA are the strict lower and upper part of A, respectively. Then the splitting of the matrix A~ is given as(31)A~=D~A~-E~A~-F~A~,with(32)D~A~=Ip⊗DA-DB⊗In,E~A~=Ip⊗EA-FBT⊗In,F~A~=Ip⊗FA-EBT⊗In.Now instead of solving the matrix equation (25), we will apply the LSMR-M algorithm to the preconditioned system(33)A~μ~-1y~=c~withy~=μ~x~,where μ~ is a preconditioner. As said, we use the SSOR preconditioner defined by(34)μ~SSOR=1ω(2-ω)(D~A~-ωE~A~)D~A~-1(D~A~-ωF~A~).We note that the np×np matrix A~ is not used explicitly. We only use the action of the linear operator A on a matrix V∈Rn×p, defined by A(V)=AV-VB. In addition, we use only matrix-by-vector products; then when using the SSOR preconditioner we have to compute, for a given V∈Rn×p, the matrix W∈Rn×p such that(35)w~=A~μ~SSOR-1v~withw~=vec⁡W,v~=vec⁡(V),or(36)w~=A~μ~SSOR-1Tv~=μ~SSORT-1A~Tv~withw~=vec⁡W,v~=vec⁡V.With setting (37)r~=μ~-1v~⟺v~=μ~r~withr~=vec⁡(R)the linear system (35) is equivalent to (38)w~=A~μ~SSOR-1v~⟺w~=A~r~.For computingR such that w~=A~μ~SSOR-1v~, we have to solve the following matrix equations:(39)(DA-ωEA)Y-Y(DB-ωFB)=ω(2-ω)V,(40)DAY-YDB=Z,(41)(DA-ωFA)R-R(DB-ωEB)=Z.The matrix equations (39) and (41) are also Sylvester matrix equations. But as was stated in [45], since the matrices involved in these equations are triangular, they are solved easily. In (39), the matrix Y can be computed from left to right and from top to bottom in each column; this corresponds to backward substitution. Equation (41) is solved in the opposite sense and this corresponds to forward substitution. Now, to compute W in (35), it is sufficient to use the action of the operator A on the matrix r~, defined by A(r~)=Ar~-r~B.To computeW in (36), first, we use the action of the operator AT on the matrix, defined by AT(v~)=ATv~-v~BT. Then, by setting(42)v~~=ATv~-v~BT,the linear system (36) is equivalent to (43)w~=μ~SSORT-1v~~⟺μ~SSORTw~=v~~.Therefore, W can be obtained by solving the following matrix equations:(44)DA-ωEATY-YDB-ωFBT=ω2-ωv~~,(45)DAY-YDB=Z,(46)(DA-ωFAT)w~-w~(DB-ωEBT)=Z.Similarly, the matrix equations (44) and (46) are also Sylvester matrix equations. But since the matrices involved in these equations are triangular, in (44), the matrix Y can be computed from right to left and from bottom to top in each row; this corresponds to forward substitution. Equation (46) is solved in the opposite sense and this corresponds to backward substitution.In Figure3, we exhibited the function log⁡10δk with(47)δk=C-(AXk-XkB)FC-(AX0-X0B)F,versus the number of iterations for LSMR-M and the SSOR-LSMR-M. Furthermore, we note that for computing the quantity RkF, (Rk is the residual matrix in kth iteration) we used the pseudocode stated in [40]. These results were obtained for ν=100,  n=300,  p=300, and ω=0.9. The initial iterative matrix was chosen as zero matrix of suitable size. As we observe by using the SSOR preconditioner the convergence rate of the LSMR-M algorithm has increased, effectively.Figure 3
Convergence history of the LSMR-M and SSOR-LSMR-M.


## 6. Conclusion

Solving the linear matrix equations is an attractive part of research. By extending the idea of LSMR method, we have proposed Algorithm2 to solve the coupled matrix equations (1) or the least-squares problem (2) over generalized symmetric matrices. By this new iterative method on the selection of special initial matrix group, we obtain the minimum Frobenius norm solutions or the minimum Frobenius norm least-squares solutions over generalized symmetric matrices. All the presented results show that the matrix LSMR iterative method is efficient to compute the solution group of the general coupled matrix equations.



---

*Source: 562529-2015-03-23.xml*