# A Note on the⊤-Stein Matrix Equation

Authors: Chun-Yueh Chiang
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

This note is concerned with the linear matrix equationX=AX⊤B+C, where the operator (·)⊤ denotes the transpose (⊤) of a matrix. The first part of this paper sets forth the necessary and sufficient conditions for the unique solvability of the solution X. The second part of this paper aims to provide a comprehensive treatment of the relationship between the theory of the generalized eigenvalue problem and the theory of the linear matrix equation. The final part of this paper starts with a brief review of numerical methods for solving the linear matrix equation. In relation to the computed methods, knowledge of the residual is discussed. An expression related to the backward error of an approximate solution is obtained; it shows that a small backward error implies a small residual. Just like the discussion of linear matrix equations, perturbation bounds for solving the linear matrix equation are also proposed in this work.


---

## Body

## 1. Introduction

Our purpose of this work is to study the so-called⊤-Stein matrix equation
(1)X=AX⊤B+C,
where A,B,C∈ℝn×n are known matrices and X∈ℝn×n is an unknown matrix to be determined. Our interest in the ⊤-Stein equation originates from the study of completely integrable mechanical systems, that is, the analysis of the ⊤-Sylvester equation
(2)AX+X⊤B=C,
where A, B, and C are matrices in ℝn×n [1, 2]. By means of the generalized inverses or QZ decomposition [3], the solvability conditions of (2) are studied in [1, 2, 4]. Suppose that the matrix pencil A-λB⊤ is regular; that is, aA+bB⊤ is invertible for some scalars a and b. The ⊤-Sylvester equation (2) can be written as
(3)(aA+bB⊤)X+X⊤(aB+bA⊤)=aC+bC⊤.Premultiplying both sides of (3) by (aA+bB⊤)-1, we have
(4)X+UX⊤V=D,
where U=(aA+bB⊤)-1, V=aB+bA⊤, and D=(aA+bB⊤)-1(aC+bC⊤). This is of the form (1). In other words, numerical approaches for solving (2) can be obtained by transforming (2) into the form of (1) and then applying numerical methods to (1) for the solution [4–6]. With this in mind, in this note we are interested in the study of ⊤-Stein matrix equation (1).Our major purpose in this work can be divided into three parts. First, we determine necessary and sufficient conditions for the unique solvability of the solution to (1). In doing so, Zhou et al. [7] transform (1) to the standard Stein equation
(5)W=AB⊤WA⊤B+AC⊤B+C
with respect to the unknown matrix W∈ℝm×n and give the following necessary condition:
(6)μν≠1,∀μ,ν∈σ(A⊤B).
Here, σ(A⊤B) is the set of all eigenvalues of A⊤B. Zhou and his coauthors show that if (5) has a unique solution, then (1) has a unique solution. However, a counterexample is provided in [7] to show that the relation (6) is only a necessary condition for the unique solvability of (1).In [4, 8], the periodic QZ (PQZ) decomposition [3] is applied to consider the necessary and sufficient conditions of the unique solvability of (1); conditions given in [8] ignore the possibility of the existence of the unique solution, while 1 is a simple root of σ(A⊤B). This condition is included in our subsequent discussion and the following remark is provided to support our observation.Remark 1.
LetA=-1 and let B=1; that is, σ(AB⊤)={-1}. It is clear that the scalar equation X=-X⊤+C has a unique solution X=C/2. But condition (6) is not satisfied by choosing μ=ν=-1.It can also be observed from Remark1 that even if (1) is uniquely solvable, it does not imply that (5) (namely, X=X+C-C) is uniquely solvable. Conditions in [4, Equation (4.6)] provided those conditions for the unique solvability of the solution to (1) via a structured algorithm. In our work, through a complete analysis of square coefficient matrices in terms of the analysis of the spectra of the matrix A⊤B, the new approach to the condition of unique solvability of the ⊤-Stein equation (1) can be obtained.Second, we present the invariant subspace method and, more generally, the deflating subspace method to solve the⊤-Stein equation. Our methods are based on the analysis of the eigeninformations for a matrix pencil. We carry out a thorough discussion to address the various eigeninformation encountered in the subspace methods. These ideas can be implemented in algorithms easily.Finally, we take full account of the error analysis of  (1). Expressions and implications such as the residual, the backward error, and perturbation bounds are derived in this work. Note that for an approximate solution Y of (1), the backward error tells us how much the matrices A, B, and C must be perturbed. An important point found in Section 5 is that a small backward error indicates a small value for the residual ℛ=Y-AY⊤B-C, but the reverse is not usually true.Beginning in Section2, we formulate the necessary and sufficient conditions for the existence of the solution of (1) directly by means of the spectrum analysis. In Section 3 we provide a deflating subspace method for computing the solution of (1). Numerical methods for solving (1) and the related residual analysis are discussed in Section 4. The associated error analysis of (1) is given in Section 5 and concluding remarks are given in Section 6.


## 2. Solvability Conditions of the Matrix Equation (1)

In order to formalize our discussion, letA⊗B be the Kronecker product of matrices A and B, let In be the n×n identity matrix, and let ∥·∥F denote the Frobenius norm.With the Kronecker product, (1) can be written as the enlarged linear system:
(7)(In2-(B⊤⊗A)𝒫)vec(X)=vec(C),
where vec(X) stacks the columns of X into a column vector and 𝒫 is the Kronecker permutation matrix [9] which maps vec(X) into vec(X⊤); that is,
(8)𝒫=∑1≤i,j≤n2‍ejei⊤⊗eiej⊤,
where ei denotes the ith column of the n2×n2 identity matrix In2. Due to the specific structure of 𝒫, it has been shown in [10, Corollary  4.3.10] that
(9)𝒫⊤(B⊤⊗A)𝒫=A⊗B⊤.
It then follows that
(10)((B⊤⊗A)𝒫)2=(B⊤⊗A)𝒫𝒫⊤(A⊗B⊤)=B⊤A⊗AB⊤,
since 𝒫2=In2 and 𝒫=𝒫⊤. Note that eigenvalues of matrices A⊤B and AB⊤ are the same. By (10) and the property of the Kronecker product [11, Theorem 4.8], we know that
(11)σ(((B⊤⊗A)𝒫)2)={λiλj∣λi,λj∈σ(A⊤B)={λ1,…,λn},1≤i,j≤n}.
That is, the eigenvalues of (B⊤⊗A)𝒫 are related to the square roots of the eigenvalues of σ(A⊤B), but from (10), no more information can be used to decide the positivity or nonnegativity of the eigenvalues of (B⊤⊗A)𝒫. A question immediately arises as to whether it is possible to obtain the explicit expression of the eigenvalues of (B⊤⊗A)𝒫, provided the eigenvalues of A⊤B are given. In the following two lemmas, we first review the periodic QZ decomposition for two matrices and then apply it to discuss the eigenvalues of (B⊤⊗A)𝒫.Lemma 2 (see [3]).
LetA and B be two matrices in ℝn×n. Then, there exist unitary matrices P,Q∈ℂn×n such that UA:=PAQ and UB:=QHB⊤PH are two upper triangular matrices.Lemma 3.
LetA and B be two matrices in ℝm×n. Then  (1)
(B⊤⊗A)𝒫=(Q⊗PH)(UA⊗UB)𝒫(QH⊗P),(2)
σ((B⊤⊗A)𝒫)={λi,±λiλj∣λi,λj∈σ(A⊤B)={λ1,…,λn},1≤i<j≤n}.Here, z denotes the principal square root of a complex number z.Proof.
Part 1 follows immediately from Lemma2 since UA=PAQ and UB=QHB⊤PH for some unitary matrices P and Q; that is,
(12)(B⊤⊗A)𝒫=(Q⊗PH)(UB⊗UA)(P⊗QH)𝒫=(Q⊗PH)(UA⊗UB)𝒫(QH⊗P).
Let the diagonal entries ofUA and UB be denoted by {aii} and {bjj}, respectively. Then, UA⊗UB is an upper triangular matrix with given diagonal entries, specified by aii and bjj. After multiplying UA⊗UB by 𝒫 from the right, the position of the entry aiibjj is changed to be in the j+n(i-1) th row and the i+n(j-1) th column of the matrix (UA⊗UB)𝒫. They are then reshuffled by a sequence of permutation matrices to form a block upper triangular matrix with diagonal entries arranged in the following order:
(13){a11b11,[0a11b22a22b110],…,[0a11bnnannb110],22a22b22,[0a22b33a33b220],…,[0annb22a22bnn0],…,22[0an-1,n-1bnnannbn-1,n-10],annbnn}.
Note that the reshuffling process is not hard to see, when n=2, UA=[a11a120a22], and UB=[b11b120b22], we have
(14)(UA⊗UB)𝒫=[a11b11a12b11a11b12a12b1200a11b22a12b220a22b110a22b12000a22b22].
However, it is conceptually simple and regular but operationally tedious to reorder (UA⊗UB)𝒫 to show this result even for n=3 and that will be left as an exercise.
By (13), it can be seen that
(15)σ((B⊤⊗A)𝒫)={aiibii,±aiiajjbiibjj,1≤i,j≤n}={λi,±λiλj,1≤i,j≤n},
where λi=aiibii∈σ(A⊤B) for 1≤i≤n.
Before demonstrating the unique solvability conditions, we need to define that a subsetΛ={λ1,…,λn} of complex numbers is said to be ⊤-reciprocal free if and only if whenever i≠j, λi≠1/λj. This definition also regards 0 and ∞ as reciprocals of each other. Then, we have the following solvability conditions of (1).Theorem 4.
The⊤-Stein matrix equation (1) is uniquely solvable if and only if the following conditions are satisfied:  (1)
the set ofσ(A⊤B)∖{-1} is ⊤-reciprocal free;(2)
-1 can be an eigenvalue of the matrix A⊤B but must be simple.Proof.
From (7), we know that the ⊤-Stein matrix equation (1) is uniquely solvable if and only if
(16)1∉σ((B⊤⊗A)𝒫).
By Lemma 3, if λ∈σ(A⊤B), then 1/λ∉σ(A⊤B). Otherwise, 1=λ·(1/λ)∈((B⊤⊗A)𝒫). On the other hand, if -1∈σ(A⊤B) and -1 is not a simple eigenvalue, then 1∈σ((B⊤⊗A)𝒫). This verifies (16) and the proof of the theorem is complete.It is worthy noting that the condition (1) of Theorem4 is contained in the condition (6) (also appear in [7, Theorem 1]), which is the necessary and sufficient conditions for the solvability of the solution to the Stein equation (5). However, as mentioned before in Remark 1 or [7, Example 2], condition (1) is just a necessary condition for unique solvability of the solution to (1). The ⊤-Stein matrix equation (1) is uniquely solvable provided that both conditions (1) and (2) of Theorem 4 are satisfying.


## 3. The Connection between Deflating Subspace and (1)

The relationship between solution of matrix equations and the matrix eigenvalue problems has been widely studied in many applications. It is famous that solution of Riccati and polynomial matrix equations can be found by computing invariant subspaces of matrices and deflating subspaces of matrix pencils [12]. This reality leads us to find some algorithms for computing solution of (1) based on the numerical computation of invariant or deflating subspaces.Given a pair ofn×n matrices A and B, recall that the function A-λB in the variable λ is said to be the matrix pencil related to the pair (A,B). For a k-dimensional subspace 𝒳∈ℂn is called a deflating subspace for the pencil A-λB if there exists a k-dimensional subspace 𝒴∈ℂn such that
(17)A𝒳⊆𝒴,B𝒳⊆𝒴;
that is,
(18)AX=YT1,BX=YT2,
where X,Y∈ℂn×k are two full rank matrices whose columns span the spaces 𝒳 and 𝒴, respectively, and matrices T1,T2∈ℂk×k. In particular, if in (18), X=Y and B=T2=I for an n×n identity matrix I, then we have the simplified formula
(19)AX=XT1.
Here, the space 𝒳 spanned by the columns of the matrix X is called an invariant subspace for A and satisfies
(20)A𝒳⊆𝒳.
One strategy to analyze the eigeninformation is to transform one matrix pencil to its simplified and equivalent form. That is, two matrix pencils A-λB and A~-λB~ are said to be equivalent if and only if there exist two nonsingular matrices P and Q such that
(21)P(A-λB)Q=A~-λB~.
In the subsequent discussion, we will use the notion ~ to describe this equivalence relation; that is, A-λB~A~-λB~.Our task in this section is to identify eigenvectors of problem (18) and then associate these eigenvectors (left and right) with the solution of (1). We begin this analysis by studying the eigeninformation of two matrices A and B, where A-λB is a regular matrix pencil.Note that for the ordinary eigenvalue problem, if the eigenvalues are different, then the eigenvectors are linearly independent. This property is also true for every regular matrix pencil and is demonstrated as follows. For a detailed proof, the reader is referred to [13, Theorem 7.3] and [14, Theorem 4.2].Theorem 5.
Given a pair ofn×n matrix A and B, if the matrix pencil A-λB is regular, then its Jordan chains corresponding to all finite and infinite eigenvalues carry the full spectral information about the matrix pencil and consist of n linearly independent vectors.Lemma 6.
LetA-λB∈ℂn×n be a regular matrix pencil. Assume that matrices Xi,Yi∈ℂn×ni, i=1,2, are full rank and satisfy the following equations:(22a)AXi=YiRi,(22b)BXi=YiSi,where Ri and Si, i=1,2, are square matrices of size ni×ni. Then  (i)
Ri-λSi∈ℂni×ni are regular matrix pencils for i=1,2,(ii)
ifσ(R1-λS1)∩σ(R2-λS2)=ϕ, then the matrix [X1X2]∈ℝn×(n1+n2) is full rank.We also need the following useful lemma.Lemma 7.
Given two regular matrix pencilsAi-λBi∈ℂni×ni, 1≤i≤2, consider the following equations with respect to U,V∈ℂn1×n2(23a)A1U=VA2,(23b)B1U=VB2.Then, if σ(A1-λB1)∩σ(A2-λB2)=ϕ, (23a) has the unique solution U=V=0.Proof.
Forn2=1, we get
(24)A1u=a2v,B1u=b2v,
where a2,b2∈ℂ, u,v∈ℂn1×1. We may without loss of generality assume that b2≠0; then A1u=(a2/b2)B1u and thus u=v=0. Now, for any n2>1, consider the generalized Schur decomposition of A2-λB2. We can assume that A2=[aij] and B2=[bij] are upper triangular matrices (i.e., aij=bij=0,1≤j<i≤n2). Denote that the ith columns of U and V are ui and vi, respectively. Thus,(25a)A1ui=∑k=1iakivk,(25b)B1ui=∑k=1ibkivk,for i=1,2,…,n2.
Ifi=1, we obtained u1=v1=0 from the above discussion. Given an integer i such that 1≤i<n2 and assume that uk=vk=0 for 1≤k≤i, we claim ui+1=vi+1=0, indeed, from (25a) and (25b), we have
(26)A1ui+1=ai+1,i+1vi+1,B1ui+1=bi+1,i+1vi+1.
Again, the case k=i+1 is following the special case n2=1. By using mathematical induction, we prove this lemma.Corollary 8.
GivenA∈ℂn×n and Λ∈ℂk×k, if σ(A)∩σ(Λ)=ϕ, then the equation with respect to U∈ℂn×k(27)AU=UΛ
has the unique solution U=0.Now we have enough tools to analyze the solution of (1) associated with some deflating spaces. We first establish an important matrix pencil; let the matrix pencil ℳ-λℒ be defined as
(28)ℳ-λℒ:=[BA⊤0-CA⊤In]-λ[In0AC⊤AB⊤]∈ℝ2n×2n;
it is clear that
(29)σ(ℳ-λℒ)=σ(BA⊤)∪σ(In-λAB⊤);
a direct calculation shows that X is a solution of the (1) if and only if
(30)ℳ[InXA⊤]=[InAX⊤]BA⊤,ℒ[InXA⊤]=[InAX⊤]
or if and only if its dual form is
(31)[-AX⊤In]ℳ=[-XA⊤In],[-AX⊤In]ℒ=AB⊤[-XA⊤In].Armed with the property given in Theorem5 and Lemma 7, we can now tackle the problem of determining how the deflating subspace is related to the solution of (1).Theorem 9.
LetA, B, and C∈ℝn×n be given in (1) and let(32a)ℳ[U1V1]=[U2V2]T1,(32b)ℒ[U1V1]=[U2V2]T2,where [UiVi] is full rank, i=1,2. Assume that the set of σ(BA⊤) is ⊤-reciprocal free. Then, one has(1)
U1=U2=0 if σ(T1-λT2)=σ(In-λAB⊤),(2)
U1 and U2 are nonsingular if T1-λT2~BA⊤-λIn. Moreover, if A is nonsingular, then X=V1U1-1A-⊤=U2-⊤V2⊤A-⊤ is the unique solution of (1).Proof.
From (32a) and (32b) we get(33a)BA⊤U1=U2T1,(33b)-CA⊤U1+V1=V2T1,(33c)U1=U2T2,(33d)AC⊤U1+AB⊤V1=V2T2.(i) It follows from (33a) and (33c) that since σ(BA⊤-λIn)∩σ(T1-λT2)=ϕ, we have U1=U2=0 by Lemma 7. (ii) It can be seen that there exist two nonsingular matrices U and V such that
(34)ℳ[0U]=[0V]T2,ℒ[0U]=[0V]T1.
Hence, together with (32a) and (32b) we have
(35)ℳ[0U1UV1]=[0U2VV2][T200T1],ℒ[0U1UV1]=[0U2VV2][T100T2].
Since σ(ℳ-λℒ)=σ(BA⊤-λIn)∪σ(In-λAB⊤) and σ(BA⊤-λIn)∩σ(In-λAB⊤)=ϕ, by Theorem 5 and Lemma 6, the matrix [0U1UV1] is nonsingular. Together with (33c), U1 and U2 are nonsingular.
LetXi=ViUi-1,i=1,2; then form (33b) and (33d)
(36)AC⊤+AB⊤X1=V2T2U1-1=X2,-CA⊤+X1=V2T1U1-1=X2BA⊤,
or
(37)AC⊤+AB⊤X1=X2,AC⊤+AB⊤X2⊤=X1⊤.
Since the set of σ(AB⊤)=σ(BA⊤) is ⊤-reciprocal free, together with
(38)X1⊤-X2-AB⊤(X1⊤-X2)⊤=0,
we get X1=X2⊤. If A is nonsingular, it is easy to verify that two matrices X1A-⊤ and X2⊤A-⊤ are both satisfying ⊤-Stein equation (1). The proof of part (ii) is complete.Remark 10.
(1) It is easily seen that[InXA⊤] and [U1V2] both span the unique deflating subspace of ℳ-λℒ corresponding to the set of σ(BA⊤). Otherwise, in part (ii) we know that T2 is nonsingular. We then are able to transform the formulae defined in (32a) and (32b) into the generalized eigenvalue problem as follows:
(39)ℳ[U1V1]=ℒ[U1V1]BA⊤.
That is, some numerical methods for the computation of the eigenspace ofℳ-λℒ corresponding to the set of σ(BA⊤) can be designed and solved (1).
(2) Since the transport of the unique solutionX of (1) is equal to the unique solution Y of the following matrix equation(40)Y=B⊤YA⊤+C⊤,
analogous to the consequences of Theorem 9, the similar results can be obtained with respect to (40) if B is nonsingular. However, we point out that (1) can be solved by computing deflating subspaces of other matrix pencils. For instance we let
(41)ℳ1-λℒ1:=[A⊤B0-C-AC⊤BIn]-λ[In00AB⊤].
Assume that the set of σ(BA⊤) is ⊤-reciprocal free; it can be shown that ℳ1[InX]=ℒ1[InX]A⊤B and it has similar results to the conclusion of Theorem 9. The unique solution X of (1) can be found by computing deflating subspaces of the matrix pencil ℳ1-λℒ1 without the assumption of the singularity of A and B.


## 4. Computational Methods for Solving (1)

Numerical methods for solving (1) have received great attention in theory and in practice and can be found in [5, 6] for the Krylov subspace methods and in [15–17] for the Smith-type iterative methods. In particular, Smith-type iterative methods are only workable in the case ρ(AB⊤)<1, where ρ(AB⊤) denotes the spectral radius of AB⊤. In the recent years, a structure algorithm has been studied for (1) [4] via PQZ decomposition, which consists of transforming into Schur form by a PQZ decomposition and then solving the resulting triangular system by way of back-substitution. In this section, we revisit these numerical methods and point out the advantages and drawbacks of all algorithms.


### 4.1. Krylov Subspace Methods

Since the⊤-Stein equation is essentially a linear system (7), we certainly can use the Krylov subspace methods to solve (7). See, for example, [5, 6] and the reference cited therein. The general idea for applying the Krylov subspace methods is defining the ⋆-Stein operator 𝒯 as 𝒯:X→X-AX⊤B and its adjoint liner operator 𝒯 as 𝒯*:Y→Y-BY⊤A such that 〈𝒯(X),Y〉=〈X,𝒯*(Y)〉. Here, X, Y∈ℝm×n and the notion 〈·,·〉 is denoted as the Frobenius inner product. Then, the iterative method based on the Krylov subspaces for (1) is as follows.(i) The conjugate gradient (CG) method [5]:
(42)Xk+1=Xk+∥Rk∥2∥Pk∥2Pk,Rk+1=C-𝒯(Xk+1)=Rk-∥Rk∥2∥Pk∥2𝒯(Pk),Dk+1=𝒯*(Rk+1)+∥Rk+1∥2∥Rk∥2Dk,
with an initial matrix X0 and the corresponding initial conditions
(43)R0=C-𝒯(X0),D0=𝒯*(R0).Note that when the solvability conditions of Theorem4 are met, the CG method is guaranteed to converge in a finite number of iterations for any initial matrix X0.


### 4.2. The Bartels-Stewart-Like Algorithm [18]

In this section, we focus on the discussion of the Bartels-Stewart algorithm, which is known to be a numerical stable algorithm, to solve⊤-Stein equations. This method is to solve (1) by means of the PQZ decomposition [18]. Its approach has been discussed in [4] and can be summarized as follows. From Lemma 3, we know that there exist two unitary matrices P and Q (see [3] for the computation procedure) such that
(44)PXQ¯-PAQ·QHX⊤P⊤·P¯BQ¯=PCQ¯.
With A^=PAQ and B^⊤=QHB⊤PH being uppertriangular, the transformed equation looks like
(45)[X^11x^12x^21x^22]-[A^11a^120a^22][X^11⊤x^21⊤x^12⊤x^22⊤][B^110b^21b^22]=[C^11c^12c^21c^22]
with X^=[X^11x^12x^21x^22]. We then have
(46)x^22-a^22x^22⊤b^22=c^22,(47)x^21-a^22x^12⊤B^11=c^21+a^22x^22⊤b^21,(48)x^12-A^11x^21⊤b^22=c^12+a^12x^22⊤b^22,(49)X^11-A^11X^11⊤B^11=C^11+a^12x^12⊤B^11+A^11x^21⊤b^21+a^12x^22⊤b^21.
Thus, the Bartels-Stewart algorithm can easily be constructed by first solving x^22 from (46), using x^22 to obtain x^12 and x^21 from (47) and (48), and then repeating the same discussion as (46)–(48) by taking advantage of the property of A^11 and B^11 being lower triangular matrices from (49).


### 4.3. Smith-Type Iterative Methods

Recently, a class method referred to as the Smith accelerative iteration has been studied for solving the Sylvester matrix equation [15]. The Smith accelerative iteration has attracted much interests in many papers (see [7, 17] and the references therein) because of its nice numerical behavior, a quadratic convergence rate, low computational cost, and high numerical reliability, despite the lack of a rigorous error analysis. Since the Sylvester matrix equation can be transformed into the Stein matrix equation with a suitable transformation, Zhou et al. proposed Smith-type iterative methods (including the Smith accelerative iteration, Smith (ℓ) iteration, and r-Smith iteration) for solving the Stein matrix equation [17] and applying Smith-type iterative methods to (1) [7]. In this section, we try to explain the Smith accelerative iteration based on the invariant subspace method and summarize the recent results from [7].Originally, the Smith-type iterative methods are developed to solve the standard Stein equation(50)X=𝒜Xℬ+𝒞,𝒜,ℬ,𝒞∈ℝn×n.
As mentioned before, the unknown X is highly related to the generalized eigenspace problems(51a)[ℬ0-𝒞I][IX]=[I00𝒜][IX]ℬ,
or
(51b)𝒜[XI][ℬ00I]=[XI][I0-𝒞𝒜].Premultiplying (51a) by the matrix [ℬ0-𝒜𝒞I] and postmultiplying (51b) by the matrix [I0-𝒞ℬ𝒜], we get
(52)[ℬ20-𝒞-𝒜𝒞ℬIn][InX]=[In00𝒜2][InX]ℬ2,𝒜2[XIn][ℬ200In]=[XIn][In0-𝒞-𝒜𝒞ℬ𝒜2].
Then, for any positive integer k>0, we obtain
(53)[ℬ2k-10-CkIn][InX]=[In00𝒜2k-1][InX]ℬ2k-1,𝒜2k-1[XIn][ℬ2k-100In]=[XIn][In0-Ck𝒜2k-1],
where the sequence {Ck} is defined by(54a)Ck=Ck-1+𝒜2k-1Ck-1ℬ2k-1,k≥1,(54b)C0=𝒞.The explicit expression of Ck is given as follows:
(55)Ck=∑i=12k-1‍𝒜i𝒞ℬi.
Under the condition ρ(𝒜)ρ(ℬ)<1, it is easy to see that {Ck} is convergent, and
(56)limsupk→∞∥X-Ck∥2k≤ρ(𝒜)ρ(ℬ);
that is, Ck converges quadratically to X as k→∞. This iterative method (54a) and (54b) is called the Smith accelerative iteration [15]. In recent years, some modified iterative methods are so-called Smith-type iteration, which are based on Smith iteration and improve its speed of convergence. See, for example, [17] and the references cited therein.Since the conditionρ(𝒜)ρ(ℬ)<1 implies that the assumptions of Theorem 4 hold, (1) is equivalent to (5). We can apply the Smith iteration to the (1) with the substitution (𝒜,ℬ,𝒞)=(AB⊤,A⊤B,C+AC⊤B). One possible drawback of the Smith-type iterative methods is that they cannot always handle the case when there exist eigenvalues λ,μ∈σ(A⊤B) such that λμ=-1 even if the unique solution X exists. Based on the solvable conditions given in this work, it is possible to develop a specific technique working on the particular case and it is a subject currently under investigation.


## 4.1. Krylov Subspace Methods

Since the⊤-Stein equation is essentially a linear system (7), we certainly can use the Krylov subspace methods to solve (7). See, for example, [5, 6] and the reference cited therein. The general idea for applying the Krylov subspace methods is defining the ⋆-Stein operator 𝒯 as 𝒯:X→X-AX⊤B and its adjoint liner operator 𝒯 as 𝒯*:Y→Y-BY⊤A such that 〈𝒯(X),Y〉=〈X,𝒯*(Y)〉. Here, X, Y∈ℝm×n and the notion 〈·,·〉 is denoted as the Frobenius inner product. Then, the iterative method based on the Krylov subspaces for (1) is as follows.(i) The conjugate gradient (CG) method [5]:
(42)Xk+1=Xk+∥Rk∥2∥Pk∥2Pk,Rk+1=C-𝒯(Xk+1)=Rk-∥Rk∥2∥Pk∥2𝒯(Pk),Dk+1=𝒯*(Rk+1)+∥Rk+1∥2∥Rk∥2Dk,
with an initial matrix X0 and the corresponding initial conditions
(43)R0=C-𝒯(X0),D0=𝒯*(R0).Note that when the solvability conditions of Theorem4 are met, the CG method is guaranteed to converge in a finite number of iterations for any initial matrix X0.


## 4.2. The Bartels-Stewart-Like Algorithm [18]

In this section, we focus on the discussion of the Bartels-Stewart algorithm, which is known to be a numerical stable algorithm, to solve⊤-Stein equations. This method is to solve (1) by means of the PQZ decomposition [18]. Its approach has been discussed in [4] and can be summarized as follows. From Lemma 3, we know that there exist two unitary matrices P and Q (see [3] for the computation procedure) such that
(44)PXQ¯-PAQ·QHX⊤P⊤·P¯BQ¯=PCQ¯.
With A^=PAQ and B^⊤=QHB⊤PH being uppertriangular, the transformed equation looks like
(45)[X^11x^12x^21x^22]-[A^11a^120a^22][X^11⊤x^21⊤x^12⊤x^22⊤][B^110b^21b^22]=[C^11c^12c^21c^22]
with X^=[X^11x^12x^21x^22]. We then have
(46)x^22-a^22x^22⊤b^22=c^22,(47)x^21-a^22x^12⊤B^11=c^21+a^22x^22⊤b^21,(48)x^12-A^11x^21⊤b^22=c^12+a^12x^22⊤b^22,(49)X^11-A^11X^11⊤B^11=C^11+a^12x^12⊤B^11+A^11x^21⊤b^21+a^12x^22⊤b^21.
Thus, the Bartels-Stewart algorithm can easily be constructed by first solving x^22 from (46), using x^22 to obtain x^12 and x^21 from (47) and (48), and then repeating the same discussion as (46)–(48) by taking advantage of the property of A^11 and B^11 being lower triangular matrices from (49).


## 4.3. Smith-Type Iterative Methods

Recently, a class method referred to as the Smith accelerative iteration has been studied for solving the Sylvester matrix equation [15]. The Smith accelerative iteration has attracted much interests in many papers (see [7, 17] and the references therein) because of its nice numerical behavior, a quadratic convergence rate, low computational cost, and high numerical reliability, despite the lack of a rigorous error analysis. Since the Sylvester matrix equation can be transformed into the Stein matrix equation with a suitable transformation, Zhou et al. proposed Smith-type iterative methods (including the Smith accelerative iteration, Smith (ℓ) iteration, and r-Smith iteration) for solving the Stein matrix equation [17] and applying Smith-type iterative methods to (1) [7]. In this section, we try to explain the Smith accelerative iteration based on the invariant subspace method and summarize the recent results from [7].Originally, the Smith-type iterative methods are developed to solve the standard Stein equation(50)X=𝒜Xℬ+𝒞,𝒜,ℬ,𝒞∈ℝn×n.
As mentioned before, the unknown X is highly related to the generalized eigenspace problems(51a)[ℬ0-𝒞I][IX]=[I00𝒜][IX]ℬ,
or
(51b)𝒜[XI][ℬ00I]=[XI][I0-𝒞𝒜].Premultiplying (51a) by the matrix [ℬ0-𝒜𝒞I] and postmultiplying (51b) by the matrix [I0-𝒞ℬ𝒜], we get
(52)[ℬ20-𝒞-𝒜𝒞ℬIn][InX]=[In00𝒜2][InX]ℬ2,𝒜2[XIn][ℬ200In]=[XIn][In0-𝒞-𝒜𝒞ℬ𝒜2].
Then, for any positive integer k>0, we obtain
(53)[ℬ2k-10-CkIn][InX]=[In00𝒜2k-1][InX]ℬ2k-1,𝒜2k-1[XIn][ℬ2k-100In]=[XIn][In0-Ck𝒜2k-1],
where the sequence {Ck} is defined by(54a)Ck=Ck-1+𝒜2k-1Ck-1ℬ2k-1,k≥1,(54b)C0=𝒞.The explicit expression of Ck is given as follows:
(55)Ck=∑i=12k-1‍𝒜i𝒞ℬi.
Under the condition ρ(𝒜)ρ(ℬ)<1, it is easy to see that {Ck} is convergent, and
(56)limsupk→∞∥X-Ck∥2k≤ρ(𝒜)ρ(ℬ);
that is, Ck converges quadratically to X as k→∞. This iterative method (54a) and (54b) is called the Smith accelerative iteration [15]. In recent years, some modified iterative methods are so-called Smith-type iteration, which are based on Smith iteration and improve its speed of convergence. See, for example, [17] and the references cited therein.Since the conditionρ(𝒜)ρ(ℬ)<1 implies that the assumptions of Theorem 4 hold, (1) is equivalent to (5). We can apply the Smith iteration to the (1) with the substitution (𝒜,ℬ,𝒞)=(AB⊤,A⊤B,C+AC⊤B). One possible drawback of the Smith-type iterative methods is that they cannot always handle the case when there exist eigenvalues λ,μ∈σ(A⊤B) such that λμ=-1 even if the unique solution X exists. Based on the solvable conditions given in this work, it is possible to develop a specific technique working on the particular case and it is a subject currently under investigation.


## 5. Error Analysis

Error analysis is a way for testing the stability of a numerical algorithm and evaluating the accuracy of an approximated solution. In the subsequent discussion, we want to consider the backward error and perturbation bounds for solving (1).As indicated in (44), matrices A^ and B^⊤ are both uppertriangular. We can then apply the error analysis for triangular linear systems in [19, Section 3.1] and [20] to obtain
(57)∥C^-(X^-A^X^⊤B^)∥F≤cm,nu(1+∥A^∥F∥B^∥F)∥X^∥F,
where cm,n is a content depending on the dimensions m and n and u is the unit roundoff. Since the PQZ decomposition is a stable process, it is true that
(58)∥C-(X-AX⊤B)∥F≤cm,n′u(1+∥A∥F∥B∥F)∥X∥F;
with a modest multiple cm,n′.Note that the inequality of the form (58) can serve as a stopping criterion for terminating iterations generated from the Krylov subspace methods [5, 6] and Smith-type iterative methods [15–17]. In what follows, we will derive the error associated with numerical algorithms, following the development in [20, 21].


### 5.1. Backward Error

Like the discussion of the ordinary Sylvester equations [20], the normwise backward error of an approximate solution Y of (1) is defined by(59)η(Y)≡min{ϵ:Y=(C+δC)+(A+δA)Y⊤(B+δB),22222∥δA∥F≤ϵα,∥δB∥F≤ϵβ,∥δC∥F≤ϵγ},
where α≡∥A∥F, β≡∥B∥F, and γ≡∥C∥F. Let ℛ≡δC+δAY⊤B+AY⊤δB+δAY⊤δB, which implies that ℛ=Y-AY⊤B-C. It can be seen that the residual ℛ satisfies
(60)∥ℛ∥F≤η(Y)(γ+∥Y∥Fαβ(2+η(Y))).
From (60), we know that a small backward error indeed implies a small relative residual ℛ. Since the coefficient matrices in (1) include nonlinearity, it appears to be an open problem to obtain the theoretical backward error with respect to the residual. Again, similar to the Sylvester equation discussed in [20, Section 16.2], the conditions under which a ⊤-Stein equation has a well-conditioned solution remain unknown.


### 5.2. Perturbation Bounds

Consider the perturbed equation(61)X+δX=(A+δA)(X+δX)⊤(B+δB)+(C+δC).
Let S(X)=X-AX⊤B be the corresponding ⊤-Stein operator. We then have S(δX)=δC+A(X+δX)⊤δB+δA(X+δX)⊤(B+δB). With the application of norm, it follows that
(62)∥δX∥F≤∥S-1∥F{∥δC∥F+∥δS∥F(∥X∥F+∥δX∥F)},
where ∥δS∥F≡∥A∥F∥δB∥F+∥δA∥F(∥B∥F+∥δB∥F). When ∥δS∥F is small enough so that 1≥∥S-1∥F·∥δS∥F, we can rearrange the previous result to be
(63)∥δX∥F∥X∥F≤∥S-1∥F1-∥S-1∥F·∥δS∥F(∥δC∥F∥X∥F+∥δS∥F).
With ∥C∥F=∥S(X)∥F≤∥S∥F·∥X∥F and the condition number κ(S)≡∥S∥F·∥S-1∥F, we arrive at the standard perturbation result
(64)∥δX∥F∥X∥F≤κ(S)1-κ(S)·∥δS∥F/∥S∥F(∥δC∥F∥C∥F+∥δS∥F∥S∥F).
Thus the relative error in X is controlled by those in A, B, and C, magnified by the condition number κ(S).On the other hand, we can also drop the high-order terms in the perturbation to obtain(65)δX-AδX⊤B=AX⊤δB+δAX⊤B+δC.
We then rewrite the system in terms of
(66)𝒬vec(δX)=[(X⊤B)⊤⊗ImIn⊗(AX⊤)Imn][vec(δA)vec(δB)vec(δC)],
where 𝒬=Imn-(B⊤⊗A)𝒫. Let ζ=max{∥δA∥F/∥A∥F,∥δB∥F/∥B∥F,∥δC∥F/∥C∥F}. It can be shown that
(67)∥δX∥F∥X∥F≤3Ψζ,
where Ψ=∥𝒬-1[α(X⊤B)⊤⊗ImβIn⊗(AX⊤)γImn]∥2/∥X∥F.A possible disadvantage of the perturbation bound (67), which ignores the consideration of the underlying structure of the problem, is that it overestimates the effect of the perturbation on the data. But this “universal” perturbation bound is accessible to any given matrices A, B, and C of (1).Unlike the perturbation bound (67), it is desirable to obtain a posteriori error bound by assuming δA=δB=0 and δC=X^-AX^⊤B-C in (61). This assumption gives rise to
(68)∥δX∥F∥X∥F≤∥P-1∥2∥R∥F∥X∥F.
It is true that while doing numerical computation, this bound given in (68) provides a simpler way for estimating the error of the solution of (1).


## 5.1. Backward Error

Like the discussion of the ordinary Sylvester equations [20], the normwise backward error of an approximate solution Y of (1) is defined by(59)η(Y)≡min{ϵ:Y=(C+δC)+(A+δA)Y⊤(B+δB),22222∥δA∥F≤ϵα,∥δB∥F≤ϵβ,∥δC∥F≤ϵγ},
where α≡∥A∥F, β≡∥B∥F, and γ≡∥C∥F. Let ℛ≡δC+δAY⊤B+AY⊤δB+δAY⊤δB, which implies that ℛ=Y-AY⊤B-C. It can be seen that the residual ℛ satisfies
(60)∥ℛ∥F≤η(Y)(γ+∥Y∥Fαβ(2+η(Y))).
From (60), we know that a small backward error indeed implies a small relative residual ℛ. Since the coefficient matrices in (1) include nonlinearity, it appears to be an open problem to obtain the theoretical backward error with respect to the residual. Again, similar to the Sylvester equation discussed in [20, Section 16.2], the conditions under which a ⊤-Stein equation has a well-conditioned solution remain unknown.


## 5.2. Perturbation Bounds

Consider the perturbed equation(61)X+δX=(A+δA)(X+δX)⊤(B+δB)+(C+δC).
Let S(X)=X-AX⊤B be the corresponding ⊤-Stein operator. We then have S(δX)=δC+A(X+δX)⊤δB+δA(X+δX)⊤(B+δB). With the application of norm, it follows that
(62)∥δX∥F≤∥S-1∥F{∥δC∥F+∥δS∥F(∥X∥F+∥δX∥F)},
where ∥δS∥F≡∥A∥F∥δB∥F+∥δA∥F(∥B∥F+∥δB∥F). When ∥δS∥F is small enough so that 1≥∥S-1∥F·∥δS∥F, we can rearrange the previous result to be
(63)∥δX∥F∥X∥F≤∥S-1∥F1-∥S-1∥F·∥δS∥F(∥δC∥F∥X∥F+∥δS∥F).
With ∥C∥F=∥S(X)∥F≤∥S∥F·∥X∥F and the condition number κ(S)≡∥S∥F·∥S-1∥F, we arrive at the standard perturbation result
(64)∥δX∥F∥X∥F≤κ(S)1-κ(S)·∥δS∥F/∥S∥F(∥δC∥F∥C∥F+∥δS∥F∥S∥F).
Thus the relative error in X is controlled by those in A, B, and C, magnified by the condition number κ(S).On the other hand, we can also drop the high-order terms in the perturbation to obtain(65)δX-AδX⊤B=AX⊤δB+δAX⊤B+δC.
We then rewrite the system in terms of
(66)𝒬vec(δX)=[(X⊤B)⊤⊗ImIn⊗(AX⊤)Imn][vec(δA)vec(δB)vec(δC)],
where 𝒬=Imn-(B⊤⊗A)𝒫. Let ζ=max{∥δA∥F/∥A∥F,∥δB∥F/∥B∥F,∥δC∥F/∥C∥F}. It can be shown that
(67)∥δX∥F∥X∥F≤3Ψζ,
where Ψ=∥𝒬-1[α(X⊤B)⊤⊗ImβIn⊗(AX⊤)γImn]∥2/∥X∥F.A possible disadvantage of the perturbation bound (67), which ignores the consideration of the underlying structure of the problem, is that it overestimates the effect of the perturbation on the data. But this “universal” perturbation bound is accessible to any given matrices A, B, and C of (1).Unlike the perturbation bound (67), it is desirable to obtain a posteriori error bound by assuming δA=δB=0 and δC=X^-AX^⊤B-C in (61). This assumption gives rise to
(68)∥δX∥F∥X∥F≤∥P-1∥2∥R∥F∥X∥F.
It is true that while doing numerical computation, this bound given in (68) provides a simpler way for estimating the error of the solution of (1).


## 6. Conclusion

In this note, we propose a novel approach to the necessary and sufficient conditions for the unique solvability of the solutionX of the ⊤-Stein equation for square coefficient matrices in terms of the analysis of the spectra σ(A⊤B). Solvability conditions have been derived and algorithms have been proposed in [4, 8] by using PQZ decomposition. On the other hand, one common procedure to solve the Stein-type equations is by means of the invariant subspace method. We believe that our discussion is the first which implements the techniques of the deflating subspace for solving ⊤-Stein matrix equation and might also give rise to the possibility of developing an advanced and effective solver in the future. Also, we obtain the theoretical residual analysis, backward error analysis, and perturbation bounds for measuring accurately the error in the computed solution of (1).



---

*Source: 824641-2013-08-07.xml*