1. Preliminaries
In the present paper, we are interested in the numerical solution of large-scale differential nonsymmetric Stein matrix equations on the time interval of the form:
where
A and
B are real, sparse and square matrices of size
and
respectively, and
E and
F are matrices of size
and
, respectively, with
. The initial condition is given in a factored form as
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
where
. The extended block Krylov subspace associated to
(see [
8,
9,
14,
15,
17]) is given by:
where
Algorithm 1 allows us to construct an orthonormal matrix
that is a basis of the block extended Krylov subspace
. Let
the restriction of the matrix
A to the block extended Krylov subspace
is given by
. Then, we have the following relations (see [
15])
where
and
.
Algorithm 1: The extended block Arnoldi algorithm (EBA) |
1. Inputs: A an matrix, V an matrix and m an integer. |
2. Compute the QR decomposition of , i.e., ; |
3. Set ; |
4. For |
(a) Set : first r columns of ; : second r columns of |
(b) ; ; |
(c) Orthogonalize i.e., |
* for |
* |
* ; |
* End for |
(d) Compute the decomposition of , i.e., ; |
5. End For. |
Throughout the paper, we use the following notations. The Frobenius inner product of the matrices X and Y is defined by , where denotes the trace of a square matrix Z. The associated norm is the Frobenius norm denoted by . The Kronecker product where . This product satisfies the properties: and .
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.
2. Low-Rank Approximate Solutions Method
In this section, we project the large differential equation by using the extended block Krylov subspaces
and
to obtain low-rank approximate solutions of Equation (
1).
We apply the extended block Arnoldi algorithm (Algorithm 1) to the pairs
and
, respectively, and to obtain orthonormal bases
and
and we have
where
We then consider approximate solutions of the large differential Stein matrix Equation (
1) that have the low-rank form
The matrix
is determined from the following Petrov–Galerkin orthogonality condition:
where
is the residual
Then, from (
4) and (
3), we obtain the low dimensional differential Stein equation
where
and
.
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.
Theorem 1. Let be the approximation obtained at step m by the extended block Arnoldi method, and the exact solution of low dimensional differential nonsymmetric Stein Equation (5). Then, the residual associated to satisfies the relationwhere is the matrix corresponding to the last rows of , and . Proof. Using the relations (
4) and (
2), we have
By using the following relations
we have
, since
, so
Now, as
exact solution of the low dimensional differential Stein equation
so
and since
and
are orthonormal, we obtain
where
,
and
. □
The approximate solution
can be given as a product of two matrices of low rank. Consider the singular value decomposition of the
matrix:
where
is the diagonal matrix of the singular values of
sorted in decreasing order. Let
and
be the
matrices of the first
l columns of
and
respectively, corresponding to the
l singular values of magnitude greater than some tolerance. We obtain the truncated singular value decomposition:
where
. Setting
, and
it follows that:
This is very important for large problems when one does not need to compute and store the approximation at each iteration.
The following result shows that the approximation is an exact solution of a perturbed differential Stein equation.
Theorem 2. Let be the approximate solution given by (2). Then, we havewhere Proof. By multiplying the (
5) left by
and right by
, we obtain
Using relationships
we find
So
where
and
. □
The next result states that the error satisfies also a differential Stein matrix equation.
Theorem 3. Let , the exact solution of (1), and be the low-rank approximate solution obtained at step m. The error satisfies the following differential Stein equation Proof. According to (
1) and (
4) we have
□
Next, we give an upper bound for the norm of the error
by using the 2-logarithmic norm defined by
. The logarithmic norm satisfies the following relation
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
where
Then, the error
can be expressed in the integral form as follows
By using
, we have
As
, so
□
3. Solving the Projected Differential Stein Matrix Equation
3.1. 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
where
and
.
The approximation
of
obtained at step
by Ros method is given by
where
and
solve the following Stein matrix equations
and
where
with
The Ros algorithm for solving the reduced differential Stein matrix Equation (
5) is summarized in Algorithm 2.
Algorithm 2: The 2-Rosenbrock method for solving the reduced NDSE (12) |
Input: . |
1. Choose h. |
2. ompute: |
3. Compute: |
4. Compute: |
5. For |
(a) Compute by solving Stein matrix Equation (14). |
(b) Compute by solving Stein matrix Equation (15). |
(c) Compute by (13). |
6. End For k. |
Output: |
3.2. 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
, let
be the approximation of
, where
is a solution of (
5). Then, the new approximation
of
obtained at step
by BDF method is defined by the implicit relation
where
is the step size,
and
are the coefficients of the BDF method as listed in the
Table 1.
The approximate
solves the following matrix equation
We assume that at each time
, the approximation
is factorized as a low-rank product
, where
and
, with
. We define
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.
Algorithm 3: The low-rank extended block Arnoldi method for DNSE |
Input: , choose a tolerance , a maximum number of iterations. |
1. For do |
2. Update , by Algorithm 1 (EBA) applied to |
3. Update , by Algorithm 1 (EBA) applied to |
4. Solve the low-dimensional problem (5) by BDF method or Ros method. |
5. If , stop. |
6. End If; |
7. End For (m); |
8. Using (7), the approximate solution is given by . |
Output: . |
4. 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 . The stopping criterion used for EBA-BDF method and EBA-Ros was or a maximum of iterations was achieved.
4.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
with homogeneous Dirichlet boundary conditions. The number of inner grid points in each direction was
and
for the operators
and
, respectively. The matrices
A and
B were obtained from the discretization of the operator
and
with the dimensions
and
, respectively. The discretization of the operator
and
yields matrices extracted from the
Lyapack package [
20] using the command
fdm_2d_matrix and denoted as A = fdm(
,’f_1(x,y)’,’f_2(x,y)’,’f(x,y)’). In this example,
10,000 and
10,000, respectively, and are named as
and
with
,
,
,
,
and
. For this experiment, we used
.
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.
4.2. Example 2
For the second set of experiments, we use the matrices
,
, and
from the Harwell Boeing Collection [
21]. The tolerance was set to
for the stop test on the residual. For the EBA-BDF and EBA-Ros methods, we used a constant timestep
,
In
Figure 2, the matrices
and
with dimensions
and
, respectively, and
. We plotted the Frobenius norms of the residuals
at final time
versus the number of extended block Arnoldi iterations for the EBA-BDF and EBA-Ros method.
In
Figure 3, the matrices
and
with dimensions
10,000 and
, respectively,
and
. We plotted the Frobenius norms of the residuals
at final time
versus the number of extended block Arnoldi iterations for the EBA-BDF and EBA-Ros method.
In
Table 3, we list the Frobenius residual norms at final time
and the corresponding CPU time for each method.
The numerical results are promising, showing the effectiveness of the proposed methods.
5. 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.
Author Contributions
Conceptualization, L.S., E.M.S. and H.T.A.; methodology, E.M.S.; software, L.S.; validation, L.S., E.M.S. and H.T.A.; formal analysis, L.S.; investigation, E.M.S.; writing—original draft preparation, L.S., E.M.S. and H.T.A.; writing—review and editing, L.S., E.M.S. and H.T.A. All authors have read and agreed to the published version of the manuscript.
Funding
This research received no external funding.
Acknowledgments
The authors should express their deep-felt thanks to the anonymous referees for their encouraging and constructive comments, which improved this paper.
Conflicts of Interest
The authors declare no conflict of interest.
References
- Van Dooren, P.M. Gramian Based Model Reduction of Large-Scale Dynamical Systems. Available online: https://citeseerx.ist.psu.edu/viewdoc/download;jsessionid=E53840A73CA1503CFAAF2F4EE3009768?doi=10.1.1.29.2500&rep=rep1&type=pdf (accessed on 29 June 2022).
- Zhou, B.; Lam, J.; Duan, G.R. On Smith-type iterative algorithms for the Stein matrix equation. Appl. Math. Lett. 2009, 22, 1038–1044. [Google Scholar] [CrossRef]
- Zhou, B.; Lam, J.; Duan, G.R. Toward solution of matrix equation X = Af(X) B + C. Linear Algebra Appl. 2011, 435, 1370–1398. [Google Scholar] [CrossRef]
- Bouhamidi, A.; Jbilou, K. A note on the numerical approximate solutions for generalized Sylvester matrix equations with applications. Appl. Math. Comput. 2008, 206, 687–694. [Google Scholar] [CrossRef]
- Calvetti, D.; Levenberg, N.; Reichel, L. Iterative methods for X − AXB = C. J. Comput. Appl. Math. 1997, 86, 73–101. [Google Scholar] [CrossRef]
- Datta, B. Numerical Methods for Linear Control Systems; Academic Press: Cambridge, MA, USA, 2004; Volume 1. [Google Scholar]
- Datta, B.N.; Datta, K. Theoretical and computational aspects of some linear algebra problems in control theory. Comput. Comb. Methods Syst. Theory 1986, 177, 201–212. [Google Scholar]
- Bentbib, A.H.; Jbilou, K.; Sadek, E.M. On some extended block Krylov based methods for large scale nonsymmetric Stein matrix equations. Mathematics 2017, 5, 21. [Google Scholar] [CrossRef]
- Agoujil, S.; Bentbib, A.H.; Jbilou, K.; Sadek, E.M. A minimal residual norm method for large-scale Sylvester matrix equations. Electron. Trans. Numer. Anal. 2014, 43, 45–59. [Google Scholar]
- Sadek, E.; Bentbib, A.H.; Sadek, L.; Talibi Alaoui, H. Global extended Krylov subspace methods for large-scale differential Sylvester matrix equations. J. Appl. Math. Comput. 2020, 62, 157–177. [Google Scholar] [CrossRef]
- Sadek, L.; Talibi Alaoui, H. Numerical methods for solving large-scale systems of differential equations. Ric. Mat. 2021, 1–18. [Google Scholar] [CrossRef]
- Sadek, L.; Talibi Alaoui, H. The extended block Arnoldi method for solving generalized differential Sylvester equations. J. Math. Model. 2020, 8, 189–206. [Google Scholar]
- Sadek, L.; Alaoui, H.T. Application of MGA and EGA algorithms on large-scale linear systems of ordinary differential equations. J. Comput. Sci. 2022, 62, 101719. [Google Scholar] [CrossRef]
- Bentbib, A.; Jbilou, K.; Sadek, E.M. On some Krylov subspace based methods for large-scale nonsymmetric algebraic Riccati problems. Comput. Math. Appl. 2015, 70, 2555–2565. [Google Scholar] [CrossRef]
- Heyouni, M. Extended Arnoldi methods for large low-rank Sylvester matrix equations. Appl. Numer. Math. 2010, 60, 1171–1182. [Google Scholar] [CrossRef]
- Hached, M.; Jbilou, K. Computational Krylov-based methods for large-scale differential Sylvester matrix problems. Numer. Linear Algebra Appl. 2018, 25, e2187. [Google Scholar] [CrossRef]
- Druskin, V.; Knizhnerman, L. Extended Krylov subspaces: Approximation of the matrix square root and related functions. SIAM J. Matrix Anal. Appl. 1998, 19, 755–771. [Google Scholar] [CrossRef]
- Rosenbrock, H.H. Some general implicit processes for the numerical solution of differential equations. Comput. J. 1963, 5, 329–330. [Google Scholar] [CrossRef]
- Butcher, J.C. Numerical Methods for Ordinary Differential Equations; John Wiley & Sons: Hoboken, NJ, USA, 2016. [Google Scholar]
- Penzl, T. LYAPACK A MATLAB Toolbox for Large Lyapunov and Riccati Equations, Model Reduction Problems, and Linear-Quadratic Optimal Control Problems. Available online: http://www.tu-chemintz.de/sfb393/lyapack (accessed on 29 June 2022).
- Davis, T. The University of Florida Sparse Matrix Collection, NA Digest, Vol. 97, No. 23. 7 June 1997. Available online: http://www.cise.ufl.edu/research/sparse/matrices (accessed on 10 June 2016).
| Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. |
© 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).