Matrix iteration algorithms for solving the generalized Lyapunov matrix equation

In this paper, we first recall some well-known results on the solvability of the generalized Lyapunov equation and rewrite this equation into the generalized Stein equation by using Cayley transformation. Then we introduce the matrix versions of biconjugate residual (BICR), biconjugate gradients stabilized (Bi-CGSTAB), and conjugate residual squared (CRS) algorithms. This study’s primary motivation is to avoid the increase of computational complexity by using the Kronecker product and vectorization operation. Finally, we offer several numerical examples to show the effectiveness of the derived algorithms.


Introduction
In this paper, we consider the generalized Lyapunov equation as follows: where A, N j ∈ R n×n (j = 1, 2, . . . , m), m n and C ∈ R n×n is symmetric, X ∈ R n×n is the symmetric solution of (1). The generalized Lyapunov equation (1) is related to several linear matrix equations displayed in Table 1. A large and growing amount of literature has considered the solution for these equations; see [1,2] and the references therein for an overview of developments and methods.
The generalized Lyapunov equation (1) often appears in the context of bilinear systems [3,4], stability analysis of linear stochastic systems [5,6], special linear stochastic differen-  tial equations [7] and other areas. For example, we discuss the origin of Eq. (1) in bilinear systems. The bilinear system is an interesting subclass of nonlinear control systems that naturally occurs in some boundary control dynamics [6]. The bilinear control system has been studied by scholars for many years and has the following state space representation: : where t is the time variable, x(t) ∈ R n , u(t) ∈ R m , y(t) ∈ R n are the stable, input and output vectors, respectively, u j (t) is the jth component of u(t). B ∈ R n×m , and A, N j , C are defined in (1).
For the bilinear control system (2), define Using the concept of reachability in [3,8], the reachability corresponding to (2) is where P is the solution of (1).
Moreover, the generalized Lyapunov equation (1) has wide applications in PDEs. Consider the heat equation subjected to mixed boundary conditions [9] x t = x in , n · ∇x = u(x -1) on 1 , where 1 , 2 , 3 and 4 are the boundaries of . For example, for a simple 2 × 2 mesh, the state vector x = [x 11 , x 21 , x 12 , x 22 ] T contains the temperatures at the inner points and the Laplacian is approximated via with meshsize h = 1/3. If the Robin condition is imposed on the whole boundary, then we have x 01 ≈ x 11hu(x 11 -1), . . .
Altogether this leads to the bilinear systeṁ where E j = e j e T j with canonical unit vector e j ∈ R 2 , and Thus, the optimal control problem of (4) reduces to the bilinear control system (2) and we ultimately need solve the generalized Lyapunov equation: Therefore, considering the important applications of the generalized Lyapunov equation (1), many researchers pay much attention to study the solution for this equation in recent years. Damm showed the direct method to solve the generalized Lyapunov equation [9]. Fan [11]. Based on the recent results, we mainly discuss the matrix iteration algorithms for the generalized Lyapunov equation (1).
The rest of the paper is organized as follows. In Sect. 2, we recall some known results on the generalized Lyapunov equation's solvability and rewrite this equation into the generalized Stein equation by using Cayley transformation. In Sect. 3, we present the matrix versions and variant forms of the BICR, Bi-CGSTAB, and CRS algorithms. In Sect. 4, we offer several numerical examples to test the effectiveness of the derived algorithms. In Sect. 5, we draw some concluding remarks.
Throughout this paper, we shall adopt the following notations. R m×n and Z + stand for the set of all m × n real matrices and positive integers. For A = (a ij ) = (a 1 , a 2 , . . . , a n ) ∈ R m×n , the symbol vec(A) is a vector defined by vec(A) = (a T 1 , a T 2 , . . . , a T n ) T . A T and A represent the transpose and 2-norm of matrix A, respectively. The symbol A ≥ 0 means that A is symmetric positive semi-definite. For B ∈ R m×n , the Kronecker product and inner product of A and B are defined by A ⊗ B = (a ij B) and A, B = tr(B T A). The open right-half and left-half planes are denoted by C + and C -, respectively.

Solvability and Cayley transformation 2.1 Solvability of the generalized Lyapunov equation
This section introduces the solvability for the generalized Lyapunov equation (1).
Denote σ (T) ∈ C by the spectrum of a linear operator T and ρ(T) = max{|λ||λ ∈ σ (T)} by the spectral radius. Define the linear matrix operators L A and : R n×n → R n×n by Obviously, (X) ≥ 0 when X ≥ 0. Therefore, using Theorem 3.9 in [6], we immediately get the generalized Lyapunov equation's stability result. Theorem 2.1 Let A ∈ R n×n and be positive. The following conclusions are equivalent: (e) σ (L A (X)) ⊂ Cand ρ(L -1 A (X) (X)) < 1, where the linear matrix operators L A and are defined by (5).

Cayley transformation for (1)
In this section, we introduce Cayley transformation for the generalized Lyapunov equation.
It is well known that Cayley transformation is a link between the classical Lyapunov and Stein equations. Fan et al. have shown that the stability of the Lyapunov and Stein equations is different. Naturally, we wonder if the stability is different and the counterparty method has other effects. It is verified in Sect. 4 that our iteration methods are more efficient after applying Cayley transformation to the generalized Lyapunov equation. We first recall the definition of Cayley transformation. Next, we show that the generalized Lyapunov equation can be changed to the generalized Setin equation after Cayley transformation. (1), take the positive parameter γ such that the matrices (γ I + A) and (γ I + A T ) are both nonsingular. Then (1) is equivalent to the generalized Stein equation

Theorem 2.2 For the generalized Lyapunov equation
Proof Introducing the positive parameter γ to (1), we get Since (γ I + A) and (γ I + A T ) are both nonsingular, premultiplying (γ I + A) -1 and postmultiplying (γ I + A T ) -1 on both sides to (7) The generalized Stein equation can be rewritten as Hence, it is not difficult to derive the following relation between A 1 and A 2 : is the preconditioning matrix and the corresponding generalized Stein equation is the preconditioning system.
By Remark 2.2, the operator P can be defined as In Sect. 3, we apply this operator to derive the variant forms of the BICR, Bi-CGSTAB, and CRS algorithms, respectively. The iteration methods are efficient. Numerical examples address this point in Sect. 4.

Iteration algorithms
This section presents the matrix versions and variant forms of the BICR, Bi-CGSTAB, and CRS algorithms in three subsections, respectively.

BICR algorithm
The BiCR method [12] has been proposed as a generalization of the conjugate residual (CR) [13] method for nonsymmetric matrices. Recently, Abe et al. designed BiCR for symmetric complex matrices (SCBiCR) and analyzed the factor in the loss of convergence speed [14]. It is easy to see that the BICR algorithm cannot be directly used to solve the generalized Lyapunov equation. Naturally, one can convert this matrix equation into the linear system through Kronecker product and vectorization operators. However, this makes the computational cost especially expensive. When the matrix order becomes larger, as the computer memory is limited, it is hard to implement in practice.
Therefore, we need to modify the BICR algorithm and ensure that the calculation cost is relatively cheap. In this subsection, we propose the matrix version of the BICR algorithm (Algorithm 1). Then we show the variant version of the BICR algorithm (Algorithm 2).
Using the iteration schemes of Algorithm 1 and Algorithm 2, we can directly solve the generalized Lyapunov equation. Further, we show the bi-orthogonal properties and convergent analysis of Algorithm 1 by Theorem 3.1 and Theorem 3.2.
For the proof of Theorem 3.1, refer to the Appendix.

Algorithm 1:
Matrix form of the BICR algorithm , , Theorem 3.2 For Algorithm 1, the relative residual error has the following property: Proof Using Theorem 3.1, we have Hence, the proof of Theorem 3.2 is completed.

Algorithm 2: The variant form of the BICR algorithm
Input: Choose initial matrices X(1) ∈ R n×n , S(1) ∈ R n×n ; Output: X with R(k) ≤ tol; (1)), k = 1; step 2: While R(k) > tol do: Remark 3.1 In terms of Theorem 3.2, the property R(k + 1) ≤ R(k) ensures that Algorithm 1 possesses fast and smooth convergence. Sonneveld [15] has shown a variant of BiCG, referred to the conjugate gradient squared (CGS). Van der Vorst [16] has derived one of the most successful variants of BiCG, known as the Bi-CGSTAB method. The Bi-CGSTAB algorithm is an effective algorithm for solving large sparse linear systems [16,17]. Chen et al. [18] proposed a flexible version of the BiCGStab algorithm for solving the linear system. It is easy to see that the Bi-CGSTAB algorithm cannot be directly used to solve the generalized Lyapunov equation. Similarly, we need to modify the Bi-CGSTAB algorithm to the matrix version. The matrix version of the Bi-CGSTAB algorithm is summarized in Algorithm 3. The variant form of the BICR algorithm is shown in Algorithm 4. Viewing the iteration schemes, we can be seen that Algorithm 3 is a simple matrix form of the Bi-CGSTAB algorithm. Hence, Algorithm 3 has the same properties as the Bi-CGSTAB algorithm. Algorithm 4 is an improved version of the Bi-CGSTAB algorithm, which has high computing efficiency. This point has been addressed by numerical examples in Sect. 4.

CRS algorithm
Zhang et al. proposed the conjugate residual squared (CRS) method in [19,20] to solve the linear system. The CRS algorithm is mainly aimed to avoid using the transpose of A in the BiCR algorithm and get faster convergence for the same computational cost [19]. Recently, Ma et al. [21] used the matrix CRS iteration method to solve a class of coupled Sylvester-transpose matrix equations. Later, they extended two mathematical equivalent X(k + 1) = X(k) + α(k + 1)P(k + 1) + ω(k + 1)S(k + 1).
forms of the CRS algorithm to solve the periodic Sylvester matrix equation by applying Kronecker product and vectorization operator [21]. In fact, in many cases, the CRS algorithm converges twice as fast as the BiCR algorithm [22,23]. The BiCR method can be derived from the preconditioned conjugate residual (CR) algorithm [24]. In exact arithmetic, they terminate after a limited number of iterations. In short, we can expect that the CRS algorithm will work well in many cases. The numerical examples in Sect. 4 are shown to address this point. It is easy to see that the CRS algorithm cannot be directly used to solve the generalized Lyapunov equation. Similarly, we need to modify the CRS algorithm to the matrix version. The matrix version of the CRS algorithm is summarized in Algorithm 5. The variant version of the CRS algorithm is shown in Algorithm 6.
Viewing the iteration schemes, it can be seen that Algorithm 5 is a simple matrix form of the CRS algorithm. Hence, Algorithm 5 has the same properties as the CRS algorithm.  The convergence result of Algorithms 2 to 6 has been summarized in Theorem 3.3. (1), if Algorithms 2 to 6 do not break down by zero division, for any initial matrix X(1) ∈ R n×n , Algorithms 2 to 6 can compute the solution of (1) within a finite number of iterations in the absence of the roundoff error.

Numerical experiments
In this section, we give several examples to show the numerical feasibility and effectiveness of Algorithm 1 (BICR), Algorithm 3 (Bi-CGSTAB algorithm), Algorithm 5 (CRS algorithm) and their improved algorithms, including Algorithm 2 (Var-BICR algorithm), Algorithm 4 (Var-Bi-CGSTAB algorithm), Algorithm 6 (Var-CRS algorithm). Set tol = 1.0e -8. The numerical behavior of iteration methods will be listed with respect to the number of iteration steps (ITs), the computing time (CPU)(s) and relative residual error (Error). All experiments are performed in Matlab (version R2017a) with double precision on a personal computer with 3.20 GHz central processing unit (Inter(R) Core(TM) i5-6500 CPU), 6.00G memory and Windows 7 operating system.

Example 4.1 Consider the generalized Lyapunov equation (1) with
Set the initial value We use Table 2 to show the error analysis for this example. Moreover, when n = 600, we use Fig. 1 to show the error analysis of Algorithms 1 to 6. By comparing with these algorithms, it is clear that the algorithms' efficiency will greatly be improved after using a Cayley transformation. The variant versions of the Bi-CGSTAB and CRS algorithms have the best efficiency.
If the Robin condition is imposed on the whole boundary, then we have where E j = e j e T j with canonical unit vector e j . The coefficient matrices N j and the columns b j of B corresponding to the left, upper, lower, and right boundaries are given by Then the above heat equation's optimal control problem reduces to solving the generalized Lyapunov equation (1). We use Table 3 to show the residual error analysis. It is obvious that the effect of the Var-Bi-CGSTAB algorithm is optimal compared with other algorithms.
Further, we use Fig. 2 to show the error analysis when n = 64. It can be seen that the variant versions of the algorithms perform better.
Since the original system is nonlinear, it is linearized by the second-order Carleman bilinear method to obtain a system of order n = m + m 2 .
The matrices A, N and b can be referred to [25]. The corresponding generalized Lyapunov equation is We use Table 4 to show the residual error analysis. Further, we use Fig. 3 to show the error analysis when n = 8. It can be seen that the Var-Bi-CGSTAB algorithm performs best.  Remark 4.1 From the three numerical examples above, it can be seen that the variant algorithms proposed in this paper will greatly improve the operating efficiency. In other words, the conjugate gradient-like methods are more efficient than the generalized Setin equation.

Conclusions
This paper has proposed the matrix versions of the BICR algorithm, Bi-CGSTAB algorithm, and CRS algorithm to solve the generalized Lyapunov equation (1). Then we have introduced the variant versions of these three algorithms. Finally, we have provided numerical examples to illustrate the feasibility and effectiveness of the derived algorithms.