Filetype PDF | Posted on 27 Jan 2023 | 3 years ago
Authentication
digital resources
268x Filetype PDF File size 0.62 MB Source: etna.math.kent.edu
Electronic Transactions on Numerical Analysis. ETNA
Volume54,pp.68–88,2021. KentStateUniversity and
Copyright © 2021, Kent State University. JohannRadonInstitute(RICAM)
ISSN1068–9613.
DOI:10.1553/etna_vol54s68
ONTHESOLUTIONOFTHENONSYMMETRICT-RICCATIEQUATION∗
† †
PETER BENNER AND DAVIDE PALITTA
Abstract. The nonsymmetric T-Riccati equation is a quadratic matrix equation where the linear part corresponds
to the so-called T-Sylvester or T-Lyapunov operator that has previously been studied in the literature. It has applications
in macroeconomics and policy dynamics. So far, it presents an unexplored problem in numerical analysis, and both
theoretical results and computational methods are lacking in the literature. In this paper we provide some sufficient
conditions for the existence and uniqueness of a nonnegative minimal solution, namely the solution with component-
wise minimal entries. Moreover, the efficient computation of such a solution is analyzed. Both the small-scale
and large-scale settings are addressed, and Newton-Kleinman-like methods are derived. The convergence of these
procedures to the minimal solution is proven, and several numerical results illustrate the computational efficiency of
the proposed methods.
Keywords. T-Riccati equation, M-matrices, minimal nonnegative solution, Newton-Kleinman method
AMSsubjectclassifications. 65F30, 15A24, 49M15, 39B42, 40C05
1. Introduction. In this paper, we consider the nonsymmetric T-Riccati operator
n×n n×n T T
R :R →R , R (X):=DX+X A−X BX+C,
T T
where A,B,C,D ∈ Rn×n,andweprovidesufficientconditions for the existence and unique-
ness of a minimal solution X ∈Rn×nto
min
(1.1) RT(X)=0.
ThesolutionofthenonsymmetricT-Riccatiequation(1.1)playsaroleinsolvingdynamics
generalized equilibrium (DSGE) problems [9, 22, 25]. “DSGE modeling is a method in
macroeconomics that attempts to explain economic phenomena, such as economic growth and
1
business cycles, and the effects of economic policy” . Equations of the form (1.1) appear in
certain procedures for solving DSGE models using perturbation-based methods [9, 25].
Taking inspiration from the (inexact) Newton-Kleinman method for standard algebraic
Riccati equations, we illustrate efficient numerical procedures for solving (1.1). Both the
small-scale and large-scale settings are addressed. In particular, in the latter framework, we
assume the matrices A and D to be such that the matrix-vector products Av and Dw require
O(n)floating point operations (flops) for any v,w ∈ Rn. This is the case, for instance, when
AandDaresparse. Moreover, we suppose B and C to be of low rank. These hypotheses
mimictheusualassumptions adopted when dealing with large-scale standard algebraic Riccati
equations; see, e.g., [2, 5, 7, 10, 13, 17, 18, 19, 23, 24, 21] and the recent survey paper [3].
The following is a synopsis of the paper. In Section 2 we present the result about the
existence and uniqueness of a minimal solution X to (1.1), namely the nonnegative solution
min
with component-wise minimal entries. A Newton-Kleinman method for the computation of
suchaX is derived in Section 3, and its convergence features are proven in Section 3.1. The
min
large-scale setting is addressed in Section 3.2, where the convergence of an inexact Newton-
Kleinman method equipped with a specific line search is illustrated. Some implementation
details of the latter procedure are discussed in Section 3.3. Several numerical results showing
∗Received July 3, 2020. Accepted September 20, 2020. Published online on November 20, 2020. Recom-
mended by Dario Bini. This work was supported by the Australian Research Council (ARC) Discovery Grant
No. DP1801038707.
†Department Computational Methods in Systems and Control Theory (CSC), Max Planck Institute for Dynamics
of Complex Technical Systems, Magdeburg, Germany ({benner,palitta}@mpi-magdeburg.mpg.de).
1https://en.wikipedia.org/wiki/Dynamic_stochastic_general_equilibrium
68
ETNA
KentStateUniversity and
JohannRadonInstitute(RICAM)
ONTHESOLUTIONOFTHENONSYMMETRICT-RICCATIEQUATION 69
the effectiveness of the proposed approaches are reported in Section 4, while our conclusions
are given in Section 5.
Throughoutthepaperweadoptthefollowingnotation: Thematrixinnerproductisdefined
as hX,Yi := trace(YTX) so that the induced norm is kXk2 = hX,Xi . I denotes the
F F n
identity matrix of order n, and the subscript is omitted whenever the dimension of I is clear
from the context. The brackets [·] are used to concatenate matrices of conforming dimensions.
In particular, a MATLAB-like notation is adopted, where [M,N] denotes the matrix obtained
by augmenting M with N. A ≥ 0 (A > 0) indicates a nonnegative (positive) matrix, that
is, a matrix whose entries are all nonnegative (positive). Clearly, A ≤ 0 (A < 0) if −A ≥ 0
(−A>0),andA≥BifA−B≥0.Moreover,werecallthatamatrixAisaZ-matrixifall
its off-diagonal entries are nonpositive. It is easy to show that a Z-matrix can be written in the
form A = sI −N, where s ∈ R and N ≥ 0. If s ≥ ρ(N), where ρ(·) denotes the spectral
radius, then A is called M-matrix.
Furthermore, we will always suppose that the following assumption holds.
ASSUMPTION 1.1. Weassumethat
• B is nonnegative (B ≥ 0) and C is nonpositive (C ≤ 0).
• I⊗D+(AT⊗I)ΠisanonsingularM-matrixwhere⊗denotestheKroneckerproduct
2 2 P P
n ×n n n
while Π ∈ R is a permutation matrix given by Π := i=1 j=1Ei,j ⊗Ej,i.
Thematrix E ∈Rn×ninAssumption1.1isthematrixwhose(i,j)-thentry is 1 while
i,j
all the others are zero.
Notice that I ⊗ D + (AT ⊗ I)Π being a nonsingular M-matrix implies that the T-
n×n n×n T
Sylvester operator S : R →R ,S (X):=DX+X A,hasanonnegativeinverse,
T T
i.e., S−1(X) ≥ 0 for X ≥ 0. For the standard Sylvester operator S : Rn×n → Rn×n,
T
S(X) := DX +XA,thisis guaranteed by assuming A, D to be nonsingular M-matrices;
see, e.g., [10, Theorem A.20]. Another important consequence of Assumption 1.1 is the
monotonicity of ST, i.e., ST(X) ≥ 0 implies X ≥ 0; see, e.g., [8].
It is not easy to analyze the impact that Assumption 1.1 has on the coefficient matrices A
andD. Nevertheless, a careful inspection of the ordering of the entries of I ⊗D+(AT ⊗I)Π
shows that if the latter is a singular M-matrix, then A ≤ 0. Indeed, every entry of A appears,
T
at least once, as an off-diagonal entry in I ⊗ D + (A ⊗I)Π, and since the off-diagonal
components of an M-matrix are nonpositive, it must hold that A ≤ 0.
2. Existence and uniqueness of a minimal solution. In this section we provide suffi-
cient conditions for the existence and uniqueness of a minimal solution X of (1.1). Our
min
result relies on the following fixed-point iteration:
X =0,
(2.1) 0
DX +XT A=XTBX −C, k≥0.
k+1 k+1 k k
THEOREM2.1. Theiterates computed by the fixed-point iteration (2.1) are such that
X ≥X , k≥0,
k+1 k
and, if there exists a nonnegative matrix Y such that R (Y ) ≥ 0, then X ≤ Y for any
T k
k ≥ 0. Moreover, under this assumption, the sequence {X } converges, and its limit is the
minimal nonnegative solution X to (1.1). k k≥0
min
ETNA
KentStateUniversity and
JohannRadonInstitute(RICAM)
70 P. BENNERANDD.PALITTA
Proof. We first show that X ≥X foranyk ≥0byinductiononk. Fork = 0,we
k+1 k
have X = S−1(−C) ≥ 0 = X asC ≤ 0. WenowassumethatX¯ ≥ X¯ for a certain
1 T 0 k k−1
¯
k > 0, and we show that X¯ ≥X¯.Wehave
k+1 k
X¯ =S−1(XTBX¯−C)
k+1 T ¯ k
k
=S−1(XTBX¯)+S−1(−C)=S−1(XTBX¯)+X +X¯−X¯
T ¯ k T T ¯ k 1 k k
k k
=S−1(XTBX¯)+X +X¯−S−1(XT BX¯ −C)
T ¯ k 1 k T ¯ k−1
k k−1
=S−1(XTBX¯−XT BX¯ )+X¯.
T ¯ k ¯ k−1 k
k k−1
Clearly, XT ≥ XT , as X¯ ≥ X¯ bytheinduction hypothesis. Therefore, recalling that
¯ ¯ k k−1
k k−1
B≥0,wehave
XTBX¯−XT BX¯ ≥0,
¯ k ¯ k−1
k k−1
so that X¯ ≥X¯.
k+1 k
Wenowsupposethatthereexistsanonnegative Y such that RT(Y) ≥ 0, and we show
that X ≤ Y for any k ≥ 0 by induction on k once again. The result is straightforward for
k
¯
k = 0 as X = 0. We now assume that X¯ ≤ Y for a certain k > 0, and we show that
0 k
X¯ ≤Y.SinceX¯ ≤Y andB ≥0,XTBX¯ ≤YTBY sothat−XTBX¯ ≥−YTBY.
k+1 k ¯ k ¯ k
Thus, we can write k k
0 ≤ DY +YTA−YTBY +C ≤DY +YTA−XTBX¯+C,
¯ k
k
and since by definition −XTBX¯ +C = −DX¯ −XT A,weget
¯ k k+1 ¯
k k+1
0 ≤ DY +YTA−DX¯ −XT A.
k+1 ¯
k+1
This means that ST(Y −X¯ ) ≥ 0, which implies Y ≥ X¯ thanks to the monotonicity
k+1 k+1
of S .
T
In conclusion, {X } is a nondecreasing, nonnegative sequence bounded from above,
k k≥0
therefore it has a finite limit lim X = X ≥ 0. Taking the limit on both sides
k→+∞ k min
of (2.1) shows that X is a solution of the equation R (X) = 0. Moreover, X is the
min T min
minimal nonnegative solution, as we have proven that X ≤Y foranynonnegative Y such
min
that R (Y) ≥ 0.
T
Asimilar result has been shown in [15, Theorem 2.3] for the (standard) nonsymmetric
Riccati equation.
3. The(inexact)Newton-Kleinmanmethod. Duetoitspossibleslowconvergencerate,
the fixed-point iteration (2.1) may not be attractive for the actual computation of the minimal
solution X , and a Newton-Kleinman-like method can be more effective for this task.
min
Thek-th iteration of the Newton method is defined as
′
R [X](X −X )=−R (X ),
T k+1 k T k
′
where R [X] denotes the Fréchet derivative of R at X. For the nonsymmetric T-Riccati
T T
operator, we have
′ T T T T T
R [X](Y)=DY +Y A−Y BX−X BY =(D−X B)Y +Y (A−BX),
T
and therefore the (k + 1)-st iterate of the Newton method is the solution of the T-Sylvester
equation
(3.1) (D−XTB)X +XT (A−BX )=−XTBX −C.
k k+1 k+1 k k k
ETNA
KentStateUniversity and
JohannRadonInstitute(RICAM)
ONTHESOLUTIONOFTHENONSYMMETRICT-RICCATIEQUATION 71
Depending on the problem size n, different state-of-the-art methods can be employed for the
solution of equations (3.1); see, e.g., [11, 12].
3
If n is moderate, say n ≤ O(10 ), then dense methods based on some decomposition
of the coefficient matrices can be employed to solve the T-Sylvester equations in (3.1). For
instance, in [11, Section 3] an algorithm based on the generalized Schur decomposition
of the pair (D,AT) is presented for efficiently solving a T-Sylvester equation of the form
DX+XTA=C.
If the problem dimension does not allow for dense matrix operations, then equations (3.1)
must be solved iteratively. The iterative solution of the T-Sylvester equations may introduce
someinexactness in the Newton scheme leading to the so-called inexact Newton-Kleinman
method and affecting the convergence features of the latter. By using tools similar to the ones
presented in [5], in Section 3.2 we show how a specific line search guarantees the convergence
of the inexact Newton method.
However, we first need to guarantee that the sequence {X } generated by (3.1) is
k k≥0
well-defined and that it converges to X ; this is the topic of the next section.
min
3.1. A convergence result. In this section we prove the convergence properties of the
Newton-Kleinmanmethod(3.1). To this end, we first recall a couple of classic results about
M-matrices; see, e.g., [8, Chapter 6].
LEMMA3.1. LetAbeaZ-matrix. ThenAisanonsingularM-matrixifandonlyifthere
exists a nonnegative vector v such that Av > 0. Moreover, if A is a nonsingular M-matrix and
B≥AisaZ-matrix,thenB isalsoanonsingularM-matrix.
Toproveconvergence of the Newton method to the minimal nonnegative solution X
min
to (1.1), we also need the following lemma.
LEMMA3.2. Assumethatthereexists a matrix Y such that RT(Y) > 0. Then it follows
that I ⊗ (D −XT B)+((A−BXmin)T ⊗I)ΠisanonsingularM-matrix.
min
Proof. Since Assumption 1.1 holds, we have that I ⊗D +(AT ⊗I)Π = rI 2 −N, with
n
N≥0,r>ρ(N),andwecanwrite
I ⊗(D−XT B)+((A−BX )T ⊗I)Π
min min
=I⊗D+(AT⊗I)Π−(I⊗XT B+(BX )T⊗I)Π)
min min
=rI−(N+(I⊗XT B+(BX )T⊗I)Π),
| min {z min }
≥0
as B, X ≥ 0. Notice that X ≥ 0 exists since the hypotheses of Theorem 2.1 are
min min
fulfilled. Therefore, I ⊗ (D − XT B)+((A−BX )T ⊗I)ΠisaZ-matrix.
min min
Moreover,
T T
(D−X B)(Y −X )+(Y −X ) (A−BX )
min min min min
=DY −XT BY −DX +XT BX
min min min min
+YTA−YTBX −XT A+XT BX .
min min min min
Since R (X ) = 0, it follows that −DX −XT A+XT BX =C. Moreover,
T min min min min min
adding and subtracting Y TBY , we get
T T
(D−X B)(Y −X )+(Y −X ) (A−BX )
min min min min
=R (Y)+(Y −X )TB(Y −X ).
T min min
To conclude, we notice that Y − X ≥0, as X is the minimal solution to (1.1) and
min min
RT(Y)>0.Therefore,
(D−XT B)(Y −X )+(Y −X )T(A−BX )≥R (Y)>0.
min min min min T
no reviews yet
Please Login to review.
