A Note on the ⊤ -Stein Matrix Equation

This note is concerned with the linear matrix equation 𝑋 = 𝐴𝑋 ⊤ 𝐵 + 𝐶 , 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 𝑋 . 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.


Introduction
Our purpose of this work is to study the so-called ⊤-Stein matrix equation where , , ∈ R × are known matrices and ∈ R × 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 where , , and are matrices in R × [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 − ⊤ is regular; that is, + ⊤ is invertible for some scalars and . The ⊤-Sylvester equation (2) can be written as Premultiplying both sides of (3) by ( + ⊤ ) −1 , we have where = ( + ⊤ ) −1 , = + ⊤ , and = ( + ⊤ ) −1 ( + ⊤ ). 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][5][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 with respect to the unknown matrix ∈ R × and give the following necessary condition: ] ̸ = 1, ∀ , ] ∈ ( ⊤ ) .
Here, ( ⊤ ) is the set of all eigenvalues of ⊤ . 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 ( ⊤ ). This condition is included in our subsequent discussion and the following remark is provided to support our observation. It can also be observed from Remark 1 that even if (1) is uniquely solvable, it does not imply that (5) (namely, = + − ) 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 ⊤ , 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 of (1), the backward error tells us how much the matrices , , and must be perturbed. An important point found in Section 5 is that a small backward error indicates a small value for the residual R = − ⊤ − , but the reverse is not usually true.
Beginning in Section 2, 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.

Solvability Conditions of the Matrix Equation (1)
In order to formalize our discussion, let ⊗ be the Kronecker product of matrices and , let be the × identity matrix, and let ‖ ⋅ ‖ denote the Frobenius norm.
With the Kronecker product, (1) can be written as the enlarged linear system: where vec( ) stacks the columns of into a column vector and P is the Kronecker permutation matrix [9] which maps vec( ) into vec( ⊤ ); that is, where denotes the th column of the 2 × 2 identity matrix 2 . Due to the specific structure of P, it has been shown in [10, Corollary 4.3.10] that It then follows that since P 2 = 2 and P = P ⊤ . Note that eigenvalues of matrices ⊤ and ⊤ are the same. By (10) and the property of the Kronecker product [11,Theorem 4.8], we know that That is, the eigenvalues of ( ⊤ ⊗ )P are related to the square roots of the eigenvalues of ( ⊤ ), but from (10), no more information can be used to decide the positivity or nonnegativity of the eigenvalues of ( ⊤ ⊗ )P. A question immediately arises as to whether it is possible to obtain the explicit expression of the eigenvalues of ( ⊤ ⊗ )P, provided the eigenvalues of ⊤ 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 ( ⊤ ⊗ )P.
Lemma 2 (see [3]). Let and be two matrices in R × . Then, there exist unitary matrices , ∈ C × such that := and := ⊤ are two upper triangular matrices.
Proof. Part 1 follows immediately from Lemma 2 since = and = ⊤ for some unitary matrices and ; that is, Abstract and Applied Analysis 3 Let the diagonal entries of and be denoted by { } and { }, respectively. Then, ⊗ is an upper triangular matrix with given diagonal entries, specified by and . After multiplying ⊗ by P from the right, the position of the entry is changed to be in the + ( −1) th row and the + ( − 1) th column of the matrix ( ⊗ )P. 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: Note that the reshuffling process is not hard to see, when = 2, = [ 11 12 0 22 ], and = [ 11 12 0 22 ], we have ] .
However, it is conceptually simple and regular but operationally tedious to reorder ( ⊗ )P to show this result even for = 3 and that will be left as an exercise.
By (13), it can be seen that where = ∈ ( ⊤ ) for 1 ≤ ≤ . Before demonstrating the unique solvability conditions, we need to define that a subset Λ = { 1 , . . . , } of complex numbers is said to be ⊤-reciprocal free if and only if whenever ̸ = , ̸ = 1/ . This definition also regards 0 and ∞ as reciprocals of each other. Then, we have the following solvability conditions of (1). (1) is uniquely solvable if and only if the following conditions are satisfied:

Theorem 4. The ⊤-Stein matrix equation
(2) −1 can be an eigenvalue of the matrix ⊤ but must be simple.
It is worthy noting that the condition (1) of Theorem 4 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.

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 of × matrices and , recall that the function − in the variable is said to be the matrix pencil related to the pair ( , ). For a -dimensional subspace X ∈ C is called a deflating subspace for the pencil that is, where , ∈ C × are two full rank matrices whose columns span the spaces X and Y, respectively, and matrices 1 , 2 ∈ C × . In particular, if in (18), = and = 2 = for an × identity matrix , then we have the simplified formula Here, the space X spanned by the columns of the matrix is called an invariant subspace for and satisfies One strategy to analyze the eigeninformation is to transform one matrix pencil to its simplified and equivalent form. That is, two matrix pencils − and̃−̃are said to be equivalent if and only if there exist two nonsingular matrices and such that In the subsequent discussion, we will use the notion ∼ to describe this equivalence relation; that is, − ∼̃−̃. 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 and , where − 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 of × matrix and , if the matrix pencil − 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 linearly independent vectors. Lemma 6. Let − ∈ C × be a regular matrix pencil. Assume that matrices , ∈ C × , = 1, 2, are full rank and satisfy the following equations: where and , = 1, 2, are square matrices of size × . Then We also need the following useful lemma.
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 M − L be defined as it is clear that a direct calculation shows that is a solution of the (1) if and only if or if and only if its dual form is Armed with the property given in Theorem 5 and Lemma 7, we can now tackle the problem of determining how the deflating subspace is related to the solution of (1). Theorem 9. Let , , and ∈ R × be given in (1) and let where [ ] is full rank, = 1, 2. Assume that the set of ( ⊤ ) is ⊤-reciprocal free. Then, one has Abstract and Applied Analysis Proof. From (32a) and (32b) we get (i) It follows from (33a) and (33c) that since ( ⊤ − ) ∩ ( 1 − 2 ) = , we have 1 = 2 = 0 by Lemma 7.
(ii) It can be seen that there exist two nonsingular matrices and such that Hence, together with (32a) and (32b) we have Since the set of ( ⊤ ) = ( ⊤ ) is ⊤-reciprocal free, together with we get 1 = ⊤ 2 . If is nonsingular, it is easy to verify that two matrices 1 −⊤ and ⊤ 2 −⊤ are both satisfying ⊤-Stein equation (1). The proof of part (ii) is complete. 1 2 ] both span the unique deflating subspace of M − L corresponding to the set of ( ⊤ ). Otherwise, in part (ii) we know that 2 is nonsingular. We then are able to transform the formulae defined in (32a) and (32b) into the generalized eigenvalue problem as follows:

Remark 10. (1) It is easily seen that [ ⊤ ] and [
That is, some numerical methods for the computation of the eigenspace of M− L corresponding to the set of ( ⊤ ) can be designed and solved (1).
(2) Since the transport of the unique solution of (1) is equal to the unique solution of the following matrix equation analogous to the consequences of Theorem 9, the similar results can be obtained with respect to (40) if is nonsingular. However, we point out that (1) can be solved by computing deflating subspaces of other matrix pencils. For instance we let Assume that the set of ( ⊤ ) is ⊤-reciprocal free; it can be shown that M 1 [ ] = L 1 [ ] ⊤ and it has similar results to the conclusion of Theorem 9. The unique solution of (1) can be found by computing deflating subspaces of the matrix pencil M 1 − L 1 without the assumption of the singularity of and . (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][16][17] for the Smithtype iterative methods. In particular, Smith-type iterative methods are only workable in the case ( ⊤ ) < 1, where ( ⊤ ) denotes the spectral radius of ⊤ . 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.

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 T as T : → − ⊤ and its adjoint liner operator T as T * : → − ⊤ such that ⟨T( ), ⟩ = ⟨ , T * ( )⟩. Here, , ∈ R × and the notion ⟨⋅, ⋅⟩ is denoted as the Frobenius inner product. Then, the iterative method based on the Krylov subspaces for (1) is as follows.

6
Abstract and Applied Analysis (i) The conjugate gradient (CG) method [5]: with an initial matrix 0 and the corresponding initial conditions Note that when the solvability conditions of Theorem 4 are met, the CG method is guaranteed to converge in a finite number of iterations for any initial matrix 0 . [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 and (see [3] for the computation procedure) such that ]. We then havê

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 Smithtype iterative methods (including the Smith accelerative iteration, Smith (ℓ) iteration, and -Smith iteration) for solving the Stein matrix equation [17] and applying Smithtype 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 As mentioned before, the unknown is highly related to the generalized eigenspace problems Then, for any positive integer > 0, we obtain where the sequence {C } is defined by The explicit expression of is given as follows: Under the condition (A) (B) < 1, it is easy to see that { } is convergent, and lim sup that is, converges quadratically to as → ∞. 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 (A) (B) < 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 (A, B, C) = ( ⊤ , ⊤ , + ⊤ ). One possible drawback of the Smith-type iterative methods is that they cannot always handle the case when there exist eigenvalues , ∈ ( ⊤ ) such that = −1 even if the unique solution 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.

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), matriceŝand̂⊤ are both uppertriangular. We can then apply the error analysis for triangular linear systems in [19, Section 3.1] and [20] to obtain − (̂−̂̂⊤̂) ≤ , u (1 +̂̂)̂, (57) where , is a content depending on the dimensions and and u is the unit roundoff. Since the PQZ decomposition is a stable process, it is true that with a modest multiple , . 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][16][17]. In what follows, we will derive the error associated with numerical algorithms, following the development in [20,21].

Abstract and Applied Analysis
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 , , and of (1).
Unlike the perturbation bound (67), it is desirable to obtain a posteriori error bound by assuming = = 0 and =̂−̂⊤ − in (61). This assumption gives rise to 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).

Conclusion
In this note, we propose a novel approach to the necessary and sufficient conditions for the unique solvability of the solution of the ⊤-Stein equation for square coefficient matrices in terms of the analysis of the spectra ( ⊤ ). 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).