On Some Extended Block Krylov Based Methods for Large Scale Nonsymmetric Stein Matrix Equations

In the present paper, we consider the large scale Stein matrix equation with a low-rank constant term AXB − X + EFT = 0. These matrix equations appear in many applications in discrete-time control problems, filtering and image restoration and others. The proposed methods are based on projection onto the extended block Krylov subspace with a Galerkin approach (GA) or with the minimization of the norm of the residual. We give some results on the residual and error norms and report some numerical experiments.


Introduction
In this paper, we are interested in the numerical solution of large scale nonsymmetric Stein matrix equations of the form: AXB − X + EF T = 0 (1) where A and B are real, sparse and square matrices of size n × n and s × s, respectively, and E and F are matrices of size n × r and s × r, respectively.Stein matrix equations play an important role in many problems in control and filtering theory for discrete-time large-scale dynamical systems, in each step of Newton's method for discrete-time algebraic Riccati equations, model reduction problems, image restoration techniques and other problems [1][2][3][4][5][6][7][8][9][10].
Direct methods for solving the matrix Equation (1), such as those proposed by Bartels-Stewart [11] and the Hessenberg-Schur [12] algorithms, are attractive if the matrices are of small size.For a general overview of numerical methods for solving the Stein matrix equation [1,2,13].
The Stein matrix Equation (1) can be formulated as an ns × ns large linear system using the Kronecker formulation: where vec(X) is the vector obtained by stacking all the columns of the matrix X, I n is the n-by-n identity matrix, and the Kronecker product of two matrices A and B is defined by A ⊗ B = [a ij B], where A = [a ij ].This product satisfies the properties (A ⊗ B)(C ⊗ D) = (AC ⊗ BD), (A ⊗ B) T = A T ⊗ B T and vec(AXB) = (B T ⊗ A)vec(X).Then, the matrix Equation (1) has a unique solution if and only if λµ = 1 for all λ ∈ σ(A) and µ ∈ σ(B), where σ(A) denotes the spectrum of the matrix A. Throughout the paper, we assume that this condition is satisfied.Moreover, if both A and B are Schur stable, i.e., σ(A) and σ(B) lie in the open unit disc, and then the solution of Equation ( 1) can be expressed as the following infinite matrix series: To solve large linear matrix equations, several Krylov subspace projection methods have been proposed (see, e.g., [1,[13][14][15][16][17][18][19][20][21][22][23][24] and the references therein).The main idea developed in these methods is to use a block Krylov subspace or an extended block Krylov subspace and then project the original large matrix equation onto these Krylov subspaces using a Galerkin condition or a minimization property of the obtained residual.Hence, we will be interested in these two procedures to get approximate solutions to the solution of the Stein matrix Equation (1).The rest of the paper is organized as follows.In the next section, we recall the extended block Krylov subspace and the extended block Arnoldi (EBA) algorithm with some properties.In Section 3, we will apply the Galerkin approach (GA) to Stein matrix equations by using the extended Krylov subspaces.In Section 4, we define the minimal residual (MR) method for Stein matrix equations by using the extended Krylov subspaces.We finally present some numerical experiments in Section 5.

The Extended Block Krylov Subspace Algorithm
In this section, we recall the EBA algorithm applied to (A, V), where V ∈ R n×r .The block Krylov subspace associated with (A, V) is defined as: The extended block Krylov subspace associated with the pair (A, V) is given as: The EBA Algorithm 1 is defined as follows [15,16,18,23]:

Algorithm 1. The Extended Block Arnoldi (EBA) Algorithm
(1) Inputs: A an n × n matrix, V an n × r matrix and m an integer Let Tm = V T m+1 AV m .Then, we have the following relations [25]: where E m = [0 2r×2(m−1)r , I 2r ] T is the matrix of the last 2r columns of the identity matrix I 2mr [23,25].
In the next section, we will define the GA for solving Stein matrix equations.

The Case: Both A and B Are Large Matrices
We consider here a nonsymmetric Stein matrix equation, where A and B are large and sparse matrices with r n and r s.We project the initial problem by using the extended block Krylov subspaces K e m (A, E) and K e m (B T , F) associated with the pairs (A, E) and (B T , F), respectively, and get orthonormal bases {V 1 , V 2 , ..., V m } and {W 1 , W 2 , ..., W m }.We then consider approximate solutions of the Stein matrix Equation (1) that have the low-rank form: where The matrix Y GA m is determined from the following Galerkin orthogonality condition: 4), we obtain the reduced Stein matrix equation: where .., 2mr and j = 1, 2, ..., 2mr, the solution Y m of the low-order Stein Equation ( 5) can be obtained by a direct method such as those described in [11].The following result on the norm of the residual R m allows us to stop the iterations without having to compute the approximation X GA m .
Theorem 1.Let X GA m be the approximation obtained at step m by the EBA algorithm.Then, the Frobenius norm of the residual R G m associated to the approximation X GA m is given by: where , and: The proof is similar to the one given at proposition 6 in [17].
In the following result, we give an upper bound for the norm of the error X − X GA m .
Theorem 2. Assume that A 2 < 1 and B 2 < 1, and let Y GA m be the exact solution of projected Stein matrix Equation ( 5) and X GA m be the approximate solution given by running m steps of the EBA algorithm.Then: Proof.The proof is similar to the one given at Theorem 2 in [27].
The approximate solution X GA m can be given as a product of two matrices of low rank.Consider the singular value decomposition of the 2mr × 2mr matrix: where Σ is the diagonal matrix of the singular values of Y MR m sorted in decreasing order.Let Y 1,l and Y 2,l be the 2mr × l matrices of the first l columns of Y 1 and Y 2 , respectively, corresponding to the l singular values of magnitude greater than some tolerance.We obtain the truncated singular value decomposition: l , it follows that: This is very important for large problems when one doesn't need to compute and store the approximation X m at each iteration.
The GA is given in Algorithm 2:

Algorithm 2. Galerkin Approach (GA) for the Stein Matrix Equations
(1) Inputs: A an n × n matrix, B an s × s matrix, E an n × r matrix and F an s × r matrix.
(3) For m = 1, 2, 3, ..., itermax (4) Compute V m , T A m , by Algorithm 1 applied to (A, E). (5) Compute W m , T B m , by Algorithm 1 applied to (B T , F). (6) Solve the low order Stein Equation ( 5) and compute R m F given by Equation (6) (7) if R m F ≤ tol, stop, (8) Using Equation ( 8), the approximate solution X GA m is given by In the next section, we consider the case where the matrix A is large while B has a moderate or a small size.

The Case: A Large and B Small
In this section, we consider the Stein matrix equation: where E is a matrix of size n × s with s << n.
In this case, we will consider approximations of the exact solution X as: where V m is the orthonormal basis obtained by applying the extended block Krylov subspace K e m (A, E).The orthogonality Galerkin condition gives: where R m is the m-th residual given by R m = AX m B − X m + E. Therefore, we obtain the projected Stein matrix equation: where T A = V T m AV m and E = V T m E. The next result gives a useful expression of the norm of the residual.

Theorem 3. Let Y GA
m the exact solution of the reduced Stein matrix Equation ( 11) and let X GA m = V m Y GA m be the approximate solution of Equation ( 9) with R m = R(X GA m ) the corresponding residual.Then: Proof.The residual is given by R , we have: This result is very important because it allows us to calculate the Frobenius norm of R m (X GA m ) without having to compute the approximate solution.
Next, we give a result showing that the error X − X m is an exact solution of a perturbed Stein matrix equation.
Theorem 4. Let X m be the approximate solution of Equation ( 9) obtained after m iterations of the EBA algorithm.Then: where Proof.Multiplying the Equation ( 11) from the left by V m , we obtain: As V m Ẽ = E, we get: where: We can now state the following result, which gives an upper bound for the norm of the error.
Theorem 5.If A 2 < 1 and B 2 < 1, then we have: Proof.By subtracting Equation ( 13) from Equation ( 9), we get: The error X m − X is the solution of the Stein matrix Equation ( 17) and can be expressed as: In the next section, we present projection methods based on extended block Krylov subspaces and MR property.

Minimal Residual Method for Large Scale Stein Matrix Equations
In this section, we present a MR method for solving large scale Stein matrix equations.A MR method for solving large scale Lyapunov matrix equation is given in [22].

The Case: Both A and B Are Large
Instead of using a Galerkin condition as we explained in the preceding section, we consider here approximate solutions X MR m = V m Y MR m W T m satisfying the following minimization property: We have the following result.
Theorem 6.The solution X MR m of the the minimization problem: where Y MR m solves the following low dimentional minimization problem: with E = V 1 R E and F = W 1 R F , the QR factorization of E and F, respectively.
Proof.We have: min One advantage of using the minimization approach is the fact that the projected problem (24) always has a solution that is not the case when one uses a GA.
The main problem is now how to solve the reduced order minimization problem (24).One possibility is the use of the preconditioned global conjugate gradient (PGCG) method.

The Preconditioned Global CG Method for Solving the Reduced Minimization Problem
In this section, we adopt the preconditioned conjugate gradient method (PCG) [28,29] to solve the reduced minimization problem (24).The normal equation associated with (24) is given by: where: Notice that L * m is the adjoint of the linear operator L m with respect to the Frobenius inner product is given by: and: We can decompose the matrices T we get the eigendecomposition: where and then the normal Equation ( 26) is now expressed as: where This expression suggests that one can use the first part as a preconditioner, that is, the matrix operator: It can be seen that the expression (29) corresponds to the normal equation of the following matrix operator: where Then, the preconditioned global CG algorithm is obtained by applying the preconditioner (30) to the normal equation associated with the matrix linear operator defined by Equation (31).This is summarized in Algorithm 3.

Algorithm 3. The Preconditioned Global Conjugate Gradient (PGCG) Algorithm.
( Notice that the use of the preconditioner P requires the solution, at each iteration, of a Stein equation.As the matrices D A and D B of these Stein matrix equations are diagonal matrices, this reduces the costs.
The MR Algorithm 4 for the Stein matrix equations is summarized as follows:

Algorithm 4. The Minimal Residual (MR) Method for Nonsymmetric Stein Matrix Equations
(1) Choose a tolerance tol > 0, a maximum number of itermax iterations (2) For m = 1, 2, 3, ..., itermax In this subsection, we apply the MR norm method to the nonsymmetric Stein Equation ( 9) in the case A large and B small.The approximate solution is given by: We have the following result, which is not difficult to prove.
Theorem 7. The solution of the minimization problem: where: The reduced minimization problem (33) can also be solved by using the preconditioned global CG method (PGCG), as we did for the problem (24).

Numerical Experiments
In this section, we present some numerical experiments of large and sparse Stein matrix equations.We compared EBA-MR and EBA-GA methods.For the GA and at each iteration m, we solved the projected Stein matrix equations by using the Bartels-Stewart algorithm [11].When solving the minimization reduced problem by the PGCG, we stopped the iterations when the relative norm of the residual was less than tol l = 10 −12 or when a maximum of kmax = 200 iterations was achieved.The algorithms were coded in Matlab 8.0 (2014).The stopping criterion used for EBA-MR and GA was R(X m ) F < 10 −7 or a maximum of m max = 100 iterations was achieved.
In all of the examples, the coefficients of the matrices E and F were random values uniformly distributed on [0, 1].
Example 1.In this first example, the matrices A and B are obtained from the centered finite difference discretization of the operators: on the unit square [0, 1] × [0, 1] with homogeneous Dirichlet boundary conditions.The number of inner grid points in each direction was n 0 and s 0 for the operators L A and L B , respectively.The matrices A and B were obtained from the discretization of the operator L A and L B with the dimensions n = n 2 0 and s = s 2 0 , respectively.The discretization of the operator L A (u) and L B (u) yields matrices extracted from the Lyapack package [30] using the command fdm_2d_matrix and denoted as A = fdm(n0,'f_1(x,y)','f_2(x,y)','f(x,y)').In this example, n = 10, 000 and s = 4900, respectively, and are named as A = fdm(n0, f 1 (x, y), f 2 (x, y), f (x, y)) and B = fdm(s0, g 1 (x, y), g 2 (x, y), g(x, y)) with f 1 (x, y) = −e xy , f 2 (x, y) = − sin(xy), f (x, y) = y 2 , g 1 (x, y) = −100e x , g 2 (x, y) = −12xy and g(x, y) = x 2 + y 2 .For this experiment, we used r = 3.
In Figure 1, we plotted the Frobenius norms of the residuals versus the number of iterations for the MR and the GAs.In Table 1, we compared the performances of the MR method and the GA.For both methods, we listed the residual norms, the maximum number of iteration and the corresponding execution time.Example 2. For the second set of experiments, we considered matrices from the University of Florida Sparse Matrix Collection [31] and from the Harwell Boeing Collection (http://math.nist.gov/MatrixMarket).
In Figure 2, we used the matrices A = pde2961 and B = fdm(s0, 100e x , 12xy, x 2 + y 2 ) with dimensions n = 2961 and s = 3600, respectively, and r = 3.In Figure 3, we used the matrices A=Themal and B=fdm(s 0 , e xy , sin(xy), x 2 − y 2 ) with dimensions n = 3456 and s = 6400, respectively, and r = 3.In Table 2, we compared the performances of the MR method and the GA.For both methods, we listed the residual norms, the maximum number of iterations and the corresponding execution time.
m and h B m represent the last 2r rows of the matrices T A m and T B m , respectively.Therefore, the normal Equation (25) can be written as:
Update V m , T