Sensitivity of the Solution to Nonsymmetric Differential Matrix Riccati Equation

: Nonsymmetric differential matrix Riccati equations arise in many problems related to science and engineering. This work is focusing on the sensitivity of the solution to perturbations in the matrix coefﬁcients and the initial condition. Two approaches of nonlocal perturbation analysis of the symmetric differential Riccati equation are extended to the nonsymmetric case. Applying the techniques of Fréchet derivatives, Lyapunov majorants and ﬁxed-point principle, two perturbation bounds are derived: the ﬁrst one is based on the integral form of the solution and the second one considers the equivalent solution to the initial value problem of the associated differential system. The ﬁrst bound is derived for the nonsymmetric differential Riccati equation in its general form. The perturbation bound based on the sensitivity analysis of the associated linear differential system is formulated for the low-dimensional approximate solution to the large-scale nonsymmetric differential Riccati equation. The two bounds exploit the existing sensitivity estimates for the matrix exponential and are alternative.


Introduction and Notations
In the present paper, we consider the nonsymmetric differential matrix Riccati equation /NDRE/Ẋ where the solution X(t) is a n × p real matrix and A ∈ R n×n , D ∈ R p×p , Q ∈ R n×p and S ∈ R p×n are the coefficient matrices and X 0 ∈ R n×p is a given initial value. We assume that the matrix is a nonsingular M-matrix, or an irreducible singular M-matrix. (Recall that a real square matrix A is said M-matrix if A = sI − B with B ≥ 0 and s ≥ r(B), where r(.) denotes the spectral radius. If s > r(B), the M-matrix A is nonsingular.) As a consequence, [1], A and D are both nonsingular M-matrices and can be decomposed as A = A 1 − A 2 and D = D 1 − D 2 , where A 2 , D 2 are positive and A 1 , D 1 are nonsingular M-matrices. The NDRE (1) can then be formulated aṡ X(t) + X(t)D 1 + A 1 X(t) = A 2 X(t) + X(t)D 2 + X(t)SX(t) + Q X(0) = X 0 .
The solution of NDRE (1) is given by the implicit formula [2] X(t) = e −tA 1 X 0 e −tD 1 + Let us now consider a nonsingular solution X * to the nonsymmetric algebraic Riccati equation In [1], it is proved that if H is assumed to be a nonsingular M-matrix, then the NDRE (1) has a global solution X(t), provided that the initial value X 0 satisfies the condition 0 ≤ X 0 ≤ X * , where for every matrices A, B ∈ R m×n , we write A ≤ B if a ij ≤ b ij for all i ∈ {1, . . . , m} and j ∈ {1, . . . , n}.
Nonsymmetric differential Riccati equations are related to linear boundary value problems arising in game and control theory, oscillation criterion problems for second order differential systems, variational calculus and theory of transport processes [3]. NDRE are an intermediate step in problems from singular perturbations and control theory when linear transformations are applied to reduce high-order systems to lower order or to partially decomposed systems. The properties of nonsymmetric differential Riccati equations determine the existence of the optimal open-loop strategies in Nash and Stackelberg control in game theory [4]. NDRE are induced via invariant embedding and interpretation formula from an "angularly shifted" transport model in the slab geometry [2]. Of mathematical interest, nonsymmetric differential Riccati equation describes the local coordinates of the restriction to a subset of the Lagrangian Grassmannian manifold. The most important results for NDRE are generalized in [5]. Fital and Guo, in [1], prove that for a suitable initial value X 0 , the initial value problem (1) has a nonnegative solution X(t), which converges to the stable equilibrium of (1). A closed formula, in terms of exponential of data matrices, for the general solution of (1), when S is invertible, is proposed in [6]. An analytical review of existing numerical methods to find the minimal nonnegative solution to lowdimensional NDRE (1) is given in [7], where the approximate low-dimensional solution to the large-scale case with low-dimensional right hand side is obtained after projecting (1) to low-dimensional differential equation by applying the extended block Arnoldi process.
The numerically computed solution contains errors, as a result of truncation of infinite series, round-off errors due to the finite precision machine arithmetic, error of stopping iteration procedures, etc. The computed perturbed solution can be represented as the exact solution of a slightly perturbed problem, simulating the effect of the errors mentioned above by equivalent perturbations in the data matrices. To estimate the actual error in the computed solution, it is important to find a bound of the error in the computed solution in terms of the perturbations in the data. Perturbation analysis of the nonsymmetric algebraic Riccati equation is given in [8], normwise, mixed and componentwise condition numbers, as well as residual bounds are proposed in [9][10][11]. To the best of our knowledge, the sensitivity and the conditioning of the NDRE are not yet analyzed.This work has two goals. First, we study the sensitivity of the solution X(t) to (1) to perturbations in the matrix coefficients A, D, Q, S, X 0 . It is done by adapting the bounds established for the symmetric differential Riccati equation, obtained by [12,13] to the nonsymmetric Riccati equation, which show good upper sensitivity. The perturbation bounds are very important for interpreting a numerically computed solution. The second objective is to apply the derived nonlocal perturbation bound to estimate the error of approximation in the solution when solving large-scale nonsymmetric differential Riccati equations by Krylov-type methods.
The paper is organized as follows: In Section 2, nonlocal sensitivity analysis of the nonsymmetric differential matrix Riccati Equation (1) is presented. An effective perturbation bound is proposed. In Section 3, an alternative perturbation bound based on the sensitivity analysis of the associated differential system is derived and then applied for estimating the error of approximation of the low-dimensional approximate solution. The bounds exploit existing sensitivity estimates of the matrix exponential. Numerical examples are presented in Section 4 to illustrate the theoretical results established in this work.
Throughout the paper, the following notations are used: R m×n denotes the space on m × n real matrices, . is the spectral norm M = [λ max (M M)] 1/2 , where λ max (N) is the maximum eigenvalue of the symmetric matrix N, . F is the Frobenius norm, A is the transpose of the matrix A ∈ R m×n , I n is the n × n unit matrix, and the symbol := stands for "equal by definition".

Nonlocal Perturbation Bound of NDRE
In this section, we will extend the approach proposed in [12] to the nonsymmetric differential Riccati Equation (1).
Let us denote by Z := (A, D, Q, S, X 0 ) the collection of matrix coefficients and by ∆Z := (∆Z, ∆D, ∆Q, ∆S, ∆X 0 ) the collection of equivalent perturbations in the data. The perturbation ∆Z ∈ ∆Z is continuous with ∆Z ≈ macheps φ(n) Z , where φ(n) is a low-order polynomial in n and macheps is the round-off unit of the machine arithmetic. If some of the matrix coefficients are not perturbed, then the corresponding perturbations are assumed to be zero. The perturbed nonsymmetric differential matrix Riccati equation, obtained from (1) by replacing the nominal values Z ∈ Z by Z + ∆Z, (∆Z ∈ ∆Z) is given by −(X(t) + ∆X(t))(D + ∆D) +(X(t) + ∆X(t))(S + ∆S)(X(t) + ∆X(t)) + Q + ∆Q, where X(t) + ∆X(t) is the solution to the perturbed nonsymmetric differential matrix Riccati Equation (3). For sufficiently small perturbations ∆Z ∈ ∆Z in the data Z ∈ Z, the solution X(t) + ∆X(t) to the perturbed Equation (3) exists and depends continuously on the elements of the perturbations ∆Z in the data Z. Let δ := [δ A ; δ D ; δ Q ; δ S ; δ X 0 ] ∈ R 5 + be the perturbation vector, where δ Z := ∆Z for ∆Z ∈ ∆Z.
Our aims in this section, are to extend the results obtained in [12] for the symmetric differential matrix Riccati equation to the nonsymmetric case (1) and to give a bound for the perturbations in the solution ∆X(t) ≤ f (δ, t) as a function of the perturbation vector δ.
From (3), taking into account (1), we can write the perturbation of the solution as defined for some matrix P ∈ R n×p , with M 2 (t, P) := −(∆A − X(t)∆S)P − P(∆D − ∆SX(t)) + P(S + ∆S)P and spectral norm From (4), we can state the following nonlocal perturbation bound: Let ∆X(t) be the perturbation in a solution X(t) to Equation (1), solved by a numerically stable algorithm in finite precision arithmetic according to the perturbation vector δ := (δ A , δ D , δ Q , δ S , δ X 0 ). Let us define the set Ω t where and Φ P (t, t 0 ) = e (t−t 0 )P is the fundamental matrix of equationη(t) = Pη(t) for some real matrix P. If δ satisfies the inclusion for Ω t given in (8), then the norm of the perturbation ∆X(t) is bounded by the nonlocal perturbation bound Proof. Premultiplying and postmultiplying the differential Equation (4) by the factors e −(t−τ)A c (τ) and e −(t−τ)D c (τ) respectively, and integrating with respect to τ from 0 to t, we obtain the equivalent integral form of the initial value problem (4) [2]: According to its definition, Φ P (t, t 0 ) satisfieṡ Equation (14) can be stated as where the operator Π is defined as The spectral norm of the operator Π(P)(t) in terms of (7), (9)-(11) is and the second order polynomial with a i (δ), i = 0, 1, 2, defined in (9) via (10), (11) is a Lyapunov majorant for the operator Π(.) such that In a similar way, for two arbitrary matrices (17) and (20) hold and for any positive number ρ such as then the operator Π(.) is a contraction on the ball M ρ : According to the fixed-point principle [14], the operator Equation (16) admits a solution ∆X(t) ∈ M ρ such that for we have which concludes the proof.
The perturbation bound formulated in Theorem 1 is a nonlocal sensitivity bound. The inclusion (12), (8) guarantees that there exists a solution X(t) + ∆X(t) of the perturbed Equation (3) for which the bound (13) holds. Let and with σ := 1 + X(t) and κ := σ + S (β + νσ 2 ). The estimate (13) becomes The assumption for H to be a nonsingular M-matrix implies that the closed loop-matrices are nonsingular M-matrices too [15]. This allows us, to facilitate the computation of the terms ν and β, to derive computable bounds for the spectral norm of the fundamental matrices based on the logarithmic norms

Sensitivity of Low-Dimensional Approximate Solutions to Large-Scale NDRE
For large NDRE with low-rank matrix Q, decomposed as Q = FG with F ∈ R n×s , G ∈ R p×s and s << n, we proposed in [7] to project the problem (1) onto extended Krylov subspace K m (A, F) and K m (D, G) applying the Extended block Arnoldi algorithm and to obtain the approximate solution where Y m solves the projected low-dimensional NDRĖ instead of the exact solution to (1). Here S m = W m SV m ∈ R 2ms×2ms , F m = V m F ∈ R 2ms×s , G m = W m G ∈ R 2ms×s and the block Hessenberg matrices T A m = V m AV m ∈ R 2ms×2ms and T D m = W m DW m ∈ R 2ms×2ms are obtained after transformation by the orthonormal matrices V m = V 1 , . . . , V m ∈ R n×2ms and W m = W 1 , . . . , W m ∈ R p×2ms composed of the orthonormal bases {V 1 , . . . , V m }, (V i ∈ R n×2s , i = 1, . . . .m) and {W 1 , . . . , W m }, (W i ∈ R p×2s , i = 1, . . . .m). The orthonormal bases are generated after applying the Extended block Arnoldi algorithm to the pairs (A, F) and (D, G); see [7] for more details.
The classical theory of Radon (see, e.g., [5]) states that any solution Y m (t) of the low-dimensional NDRE (24) is locally equivalent to a solution of the initial value probleṁ The solution of (25) is If the matrix Y 1,m (t) is nonsingular, the solution Y m (t) of the projected low-order nonsymmetric differential Riccati Equation (24) is represented as [16] To estimate the sensitivity of the problem (24), we consider the strategies proposed in [13]. We represent the calculated perturbed solution to (24) with collection of data , reflects the effect of round-off errors and approximation errors in the computed solution to (24).
The perturbed projected low-dimensional NDRE is The equivalent to (28) initial value problem is The perturbations ∆Ψ 11 (t), ∆Ψ 12 (t), ∆Ψ 21 (t), ∆Ψ 22 (t) are analytical functions of the data perturbations ∆Z m ∈ ∆Z m and reflect the errors in the solution Ψ(t) + ∆Ψ(t) to the perturbed linear differential system According to (26) and (30), the perturbation ∆Ψ(t) is If the matrix with ∆Y 1,m (t) and ∆Y 2,m (t) given in (31) to the perturbed projected low-dimensional Equation (28) exists. The perturbation bound of the solution Y m (t) to the NDRE (24) consists of finding an interval T : [0, t * ) such that for each t ∈ T, the matrix Y 1,m (t) + ∆Y 1,m (t) is invertible and then the perturbed solution Y m (t) + ∆Y m (t) given by (32) exists, as well as to derive a normwise bound in terms of spectral norm for the error ∆Y m (t) in the solution Y m (t) + ∆Y m (t) (32) as a function of the equivalent perturbations ∆Z m in the data coefficients Z m .
We formulate the following perturbation bound of the solution to the projected lowdimensional NDRE (24).
In order to represent the norm ∆Ψ(t) of the perturbation ∆Ψ(t) in terms of ω 1 and ω 2 by the norms of the perturbations in the data matrices T A m , T D m , F m , G m , S m and Y 0 , we consider the perturbed differential Equation (29). We have Then, taking the spectral norm, we have where δ H m := ∆H m < δ m and g(t) is an upper bound for e H m (t) , i.e., e H m (t) ≤ g(t). Some bounds for the matrix exponential e H m (t) based on Jordan and Schur matrix decompositions, logarithmic norm and power series are summarized in [18]: with constants c 0 , , and p, listed in Table 1.
Here µ(H m (t)) = λ max H m (t) + H m (t) )/2 ; ς ≥ 1 is the dimension of the maximum block in the Jordan canonical form J = Y −1 H m (t)Y of H m (t), where the matrix Y is chosen so as the condition number cond(Y) = Y Y −1 to be minimized; d ς = cos π ς+1 ; α(H m (t)) is the spectral abscissa of H m (t), i.e., the maximum real part of the eigenvalues of H m (t); T = U H H m (t)U = Λ + N is the Schur decomposition of H m (t) with unitary matrix U, chosen so as = N to be minimized, diagonal matrix Λ and N -strictly upper triangular matrix with index of nilpotency l = min{s : N s = 0}. Table 1. Values of the constants in the matrix exponential bounds.

Jordan (1) Jordan (2) Schur Log Norm Power Series
The results stated in Theorem 2 can be used to formulate a perturbation bound for the solution X(t) to the large-scale NDRE (1). (1) for which the constant matrix coefficient Q is low-rank and can decomposed as Q = FG with F ∈ R n×s , G ∈ R p×s , (s << n) be projected onto a pair of extended Krylov subspaces K m (A, F) and K m (D, G). Let X m (t) = V m Y m (t)W m ∈ R n×p be its approximate solution as stated in (23), obtained by applying the Extended block Arnoldi algorithm to the projected low-dimensional NDRĖ

Theorem 3. Let a large-scale NDRE
with solution Y m (t) and perturbation bound ∆Y m (t) ≤ f m (t, δ m ) (35) as defined in Theorem 2. A perturbation bound in terms of spectral norm for the approximate solution X(t) to (1) is given by for t ∈ [0, t * ), with t * as defined in (34) and Y 1,m (t) stated in (26), ω 1 (t, δ m ), ω(t, δ m ) stated in (33), Proof. The proof follows directly from the definition (23) of the approximate solution Y m (t) and the preservation of the spectral norm by unitary matrices.
Next, we apply the preceding results to estimate the approximation error E m = X(t) − X m (t) of the approximate solution X m (t) to the projected low-dimensional NDRE versus the exact solution to (1) from Theorem (4), which was already established by the authors in [19].

) and X(t) is an exact solution X(t) of (1).
Let us rewrite the NDRE associated with the error E m (t) where As the unitary matrices V m and W m have unit spectral norm, we have ∆ A m = T A m+1,n and ∆ D m = T D m+1,n . We notice that the fact that the term T m+1,n tends normwisely to 0 as m increases implies that the spectral norm of ∆ A m , ∆ D m decreases towards 0. This allows us to consider ∆ A m , ∆ D m as equivalent data perturbations, and Equation (43) where a 0 (δ) = ν S ; Proof. Comparing M 1,m (t, E m (t)), M 2,m (t, E m (t)) (45) to M 1 (t, P), M 2 (t, P) (6), for the coefficients a i (δ) (9) we obtain Replacing the coefficients a i (δ), i = 0, 1, 2 in the inequality of the set Ω t,m (8) and the bound f (δ, t) (13) we obtain 2 a 0 (δ)a 2 (δ) ≤ 1, or a 1 (δ)a 2 (δ) ≤ 0.25.

Numerical Examples
To illustrate the effectiveness of the bounds proposed in Theorems 1 and 2, we consider nonsymmetric differential matrix Riccati equations of type (1) on a time interval T = [0, 1], for different matrix coefficients and for several sizes. The experimental tests are performed with Matlab R2020a on an Intel processor laptop equipped with 16GB of RAM. The reference solutions X(t) to the NDRE (1) and X(t) + ∆X(t) to the perturbed NDRE (3) are computed by the backward differential formula -BDF1-Newton method, see [19] for more details. (1), constructed according to the rules given in [15]. This scheme is used in [10] to analyze the effectiveness of mixed and componentwise condition numbers, and in [11], to illustrate the validity of a condition number and backward errors of nonsymmetric algebraic Riccati equation.
The same experimental statement is used to test the accuracy of the estimate (34) and (35) from Theorem 2. The results obtained for j = 12, 10, 8, 6, 4, 2 are listed in Table 4 for α = 0 and in Table 5 for α = 5. Comparing the results for the bound (34) and (35) to these for the bound (8)- (13) given in Tables 2 and 3, it is seen that the two bounds are of the same size of the domain of validity. The bound (8)-(13) from Theorem 1 is superior to the bound (34) and (35) from Theorem 2 with respect of closeness to the estimated quantity. However, the bound (34) and (35) from Theorem 2 has the advantage that it is not related with the solution of the NDRE and hence with problems of divergence of the numerical procedure.    Table 5. Relative perturbation ρ 1 (t) = ∆X(t) X(t) and bound ρ 2 (t) = f (t) X(t) .

Experiment 2.
The size n of the matrices A, D, S, Q ∈ R n×n varies from 5 to 50. The perturbations are randomly generated following the scheme ∆Z = (rand(size(Z))/ Z ) * 10 −j for Z = A, D, S, Q and j = 12. The average values over 30 trials for the relative perturbation (8) It appears in Figure 1 that for α = 0 and j from 10 to 4, the bound ρ 2 (t) remains of the order of the estimated value ρ 1 (t).
The average values over 30 trials for the relative perturbation ρ 1 (t) = ∆X(t) X(t and the (8)
The results, obtained for the relative perturbation ρ 1 (t) = ∆X(t) X(t) and the nonlocal bound ρ 2 (t) = f (t) X(t) , f (t) = f (δ, t), (8)-(13) for t = 0.4, 0.8, 1.2, 1.6 and 2 are shown in Table 6. Table 6. Relative perturbations ρ 1 = ∆X(t) X(t) and ρ 2 = f (t) X(t) . As is seen, over all the interval of integration T ∈ 0, 2 , the perturbation bound f (δ, t) from (8)-(13) is valid, i.e., the condition (8) for existence of the bound f (δ, t) is not deteriorated. The bound f (δ, t) is a quite sharp upper bound-remains in the order of the estimated value. The perturbation bound f (δ, t) formulated in Theorem 1 is effective and could be used to estimate the sensitivity even of a large-scale NDRE.

Conclusions
In this paper, a nonlocal sensitivity analysis of the nonsymmetric differential matrix Riccati equation is presented. Two computable perturbation bounds are derived using the techniques of Fréchet derivatives, Lyapunov majorants and fixed-point principles, developed in [14]. The first bound is based on the integral form of the solution. The second one exploits the statement of the classical Radon's theory of local equivalence of the solution to the differential matrix Riccati equation to the solution of the initial value problem of the associated differential system. It has the advantage of not being related with the solution of the NDRE and hence with problems of divergence of the numerical procedure. Numerical examples show that the estimates proposed are fairly sharp for both low-dimensional and large-scale NDRE. The perturbation bound is a crucial issue of the process of numerical solution of an equation as well as a tool to evaluate the stability of the computation process. The tight perturbation bounds, proposed in the paper, allow estimation of the accuracy of the solution to a numerically solved nonsymmetric differential matrix Riccati equation. Funding: This work was partially supported by grant BG PLANTNET "Establishment of national information network genbank-Plant genetic resources".