On Some Numerical Methods for Solving Large Differential Nonsymmetric Stein Matrix Equations

: In this paper, we propose a new numerical method based on the extended block Arnoldi algorithm for solving large-scale differential nonsymmetric Stein matrix equations with low-rank right-hand sides. This algorithm is based on projecting the initial problem on the extended block Krylov subspace to obtain a low-dimensional differential Stein matrix equation. The obtained reduced-order problem is solved by the backward differentiation formula (BDF) method or the Rosenbrock (Ros) method, the obtained solution is used to build the low-rank approximate solution of the original problem. We give some theoretical results and report some numerical experiments.


Preliminaries
In the present paper, we are interested in the numerical solution of large-scale differential nonsymmetric Stein matrix equations on the time interval [t 0 , T f ] of the form: where A and B are real, sparse and square matrices of size n × n and p × p, respectively, and E and F are matrices of size n × r and p × r, respectively, with r n, p. The initial condition is given in a factored form as X 0 = Z 0Z T 0 ∈ R n×p and the matrices A and B are assumed to be large and sparse. Differential Stein matrix equations play an important role in many problems in control and filtering theory for discrete-time large-scale dynamical systems and other problems; see [1][2][3][4][5][6][7][8]. For solving the large matrix equations during the last years, several projection methods based on Krylov subspaces have been proposed, for example Sylvester matrix equations, differential Sylvester matrix equations, Riccati and others, see, e.g., [8][9][10][11][12][13][14][15][16]. The main idea employed in these methods is to project the initial problem by using extended Krylov subspaces and then apply the Petrov-Galerkin orthogonality condition.
In this work, we present an extended block Krylov method for solving large differential Stein matrix equations. Our method uses the Petrov-Galerkin projection method and the extended block Arnoldi algorithm. The problem (1) is projected onto small extended Krylov subspaces to obtain low-order differential Stein equations that are solved by a BDF method or Ros method. The approximate solution is then given as a product of matrices with low rank.
We recall the extended block Arnoldi algorithm applied to (A, V) where V ∈ R n×r . The extended block Krylov subspace associated to (A, V) (see [8,9,14,15,17]) is given by: Algorithm 1 allows us to construct an orthonormal matrix V m = [V 1 , V 2 , . . . , V m ] that is a basis of the block extended Krylov subspace K e m (A, V). Let T m,A the restriction of the matrix A to the block extended Krylov subspace K e m (A, V) is given by T m,A = V T m A V m . Then, we have the following relations (see [15]) Algorithm 1: The extended block Arnoldi algorithm (EBA)

1.
Inputs: A an n × n matrix, V an n × r matrix and m an integer.
Throughout the paper, we use the following notations. The Frobenius inner product of the matrices X and Y is defined by X, Y F = tr(X T Y), where tr(Z) denotes the trace of a square matrix Z. The associated norm is the Frobenius norm denoted by F . The Kronecker product A ⊗ B = [a i,j B]∈ R np×np where A = [a i,j ]. This product satisfies the properties: The rest of the paper is organized as follows. In Section 2, we derive the low-rank approximate solutions method and give some theoretical results. In Section 3, we apply the backward differentiation formula method and Rosenbrock method for solving reduced order problem. Section 4 is devoted to some numerical examples.

Low-Rank Approximate Solutions Method
In this section, we project the large differential equation by using the extended block Krylov subspaces K e m (A, E) and K e m (B T , F) to obtain low-rank approximate solutions of Equation (1).
We apply the extended block Arnoldi algorithm (Algorithm 1) to the pairs (A, E) and (B T , F), respectively, and to obtain orthonormal bases We then consider approximate solutions of the large differential Stein matrix Equation (1) that have the low-rank form The matrix Y m (t) is determined from the following Petrov-Galerkin orthogonality condition: where R m (t) is the residual Then, from (4) and (3), we obtain the low dimensional differential Stein equation where E m = V T m E and F m = W T m F. In the following result, we derive an expression for computation of the norm of the residual R, without having to compute matrix products with the large matrices A and B. This result allows us to reduce the cost in the proposed method when checking if the residual norm is less than some fixed tolerance.
m be the approximation obtained at step m by the extended block Arnoldi method, and Y m (t) the exact solution of low dimensional differential nonsymmetric Stein Equation (5). Then, the residual R m (t) associated to X m (t) satisfies the relation where Y m,1 (t) is the 2r × 2mr matrix corresponding to the last 2r rows of Y m (t), Y m, Proof. Using the relations (4) and (2), we have Now, as Y m (t) exact solution of the low dimensional differential Stein equatioṅ and since V m+1 and W m+1 are orthonormal, we obtain The approximate solution X m (t) 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 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: This is very important for large problems when one does not need to compute and store the approximation X m at each iteration.
The following result shows that the approximation X m (t) is an exact solution of a perturbed differential Stein equation.
Theorem 2. Let X m (t) be the approximate solution given by (2). Then, we havė where Proof. By multiplying the (5) left by V m and right by W m T , we obtain Using relationships we finḋ The next result states that the error E m (t) = X(t) − X m (t) satisfies also a differential Stein matrix equation.
Theorem 3. Let X(t), the exact solution of (1), and X m (t) be the low-rank approximate solution obtained at step m. The error E m (t) satisfies the following differential Stein equatioṅ Proof. According to (1) and (4) we havė Next, we give an upper bound for the norm of the error X − X m by using the 2-logarithmic norm defined by µ 2 (A) = 1 2 λ max (A + A T ). The logarithmic norm satisfies the following relation e tA 2 ≤ e µ 2 (A)t ; t ≥ 0.
Theorem 4. At step m of the extended block Arnoldi process, we have the following upper bound for the norm of the error, Proof. We notice that the differential Stein Equation (9) is equivalent to the linear ordinary differential equation ė rr m (t) = Aerr m (t) − b m (t), err 0 = vec(E m (t 0 )), (11) where Then, the error err m (t) can be expressed in the integral form as follows By using e tA 2 ≤ e µ 2 (A)t , we have

Rosenbrock Method
In this section, we are applying the Ros method [18] for solving the projected differential Stein matrix equation. We can write the low dimensional nonsymmetric differential Stein Equation (5) in the following form The approximation Y m,k+1 of Y m (t k+1 ) obtained at step k + 1 by Ros method is given by where P 1 and P 2 solve the following Stein matrix equations and where The Ros algorithm for solving the reduced differential Stein matrix Equation (5) is summarized in Algorithm 2.

6.
End For k.
Output: Y T f .

Backward Differentiation Formula Method
In this section, we present a BDF method [19] for solving, at each step m of the extended block Arnoldi process, the low dimensional differential Stein matrix Equation (5).
At each time t k , let Y m,k be the approximation of Y m (t k ), where Y m is a solution of (5). Then, the new approximation Y m,k+1 of Y m (t k+1 ) obtained at step k + 1 by BDF method is defined by the implicit relation where h = t k+1 − t k is the step size, α i and β are the coefficients of the BDF method as listed in the Table 1.
The approximate Y m,k+1 solves the following matrix equation We assume that at each time t k , the approximation Y m,k is factorized as a low-rank product Y m,k ≈ U m,k V T m,k , where U m,k ∈ R n×m k and V m,k ∈ R p×m k , with m k n, p.
Then, we obtain the following matrix Stein equation: The low-rank approximate solutions method by extended block Arnoldi algorithm for the large differential nonsymmetric stein matrix equations is summarized as follows in Algorithm 3.
Update V m , T m,A , by Algorithm 1 (EBA) applied to (A, E) 3.
Update W m , T m,B , by Algorithm 1 (EBA) applied to (B T , F) 4.
Solve the low-dimensional problem (5) by BDF method or Ros method.
Using (7), the approximate solution X m is given by X m ≈ Z 1,m Z T 2,m . Output: X m .

Numerical Experiments
In this section, we give some numerical examples of large nonsymmetric differential Stein matrix equations. All the experiments were performed on a computer of Intel Core i5 at 1.3 GHz and 4 GB of RAM. The Algorithm 3 were coded in Matlab2014. In all of the examples, the coefficients of the matrices E and F were random values uniformly distributed on [0, 1]. The stopping criterion used for EBA-BDF method and EBA-Ros was R(X m ) F < 10 −10 or a maximum of m max = 40 iterations was achieved.

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 p 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 p = p 2 0 , respectively. The discretization of the operator L A (u) and L B (u) yields matrices extracted from the Lyapack package [20] using the command fdm_2d_matrix and denoted as A = fdm(n 0 ,'f_1(x,y)','f_2(x,y)','f(x,y)'). In this example, n = 10,000 and p = 10,000, respectively, and are named as A = fdm(n 0 , f 1 (x, y), f 2 (x, y), f (x, y)) and B = fdm(p0, 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 = 2.
In Figure 1, we plotted the Frobenius norms of the residuals versus the number of iterations for the EBA-BDF and the EBA-Ros method. In Table 2, we compared the performances of the EBA-BDF method and the EBA-Ros. 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 use the matrices add32, pde2961, and thermal from the Harwell Boeing Collection [21]. The tolerance was set to 10 −8 for the stop test on the residual. For the EBA-BDF and EBA-Ros methods, we used a constant timestep h = 0.1, In Figure 2, the matrices A = add32 and B = pde2961 with dimensions n = 4960 and p = 2961, respectively, and r = 2. We plotted the Frobenius norms of the residuals R m (T f ) F at final time T f versus the number of extended block Arnoldi iterations for the EBA-BDF and EBA-Ros method.

Conclusions
We presented, in this paper, a new iterative method for computing numerical solutions for large-scale differential Stein matrix equations with low-rank right-hand sides. This approach is based on projection onto extended block Krylov subspaces with a Galerkin method. The approximate solutions are given as products of two low-rank matrices and allow for saving memory for large problems. The numerical experiments show that the proposed extended block Krylov-based method is effective for large and sparse matrices.