1 On the numerical solution of large-scale sparse discrete-time Riccati equations

By Heike Faßbender and Peter Benner.

Inspired by a large-scale sparse discrete-time Riccati equation which arises in a spectral factorization problem the efficient numerical solution of such Riccati equations is studied in this work. Spectral factorization is a crucial step in the solution of linear quadratic estimation and control problems. A variety of methods has been developed over the years for the computation of canonical spectral factors for processes with rational spectral densities, see, e.g., the survey [?]. One approach involves the spectral factorization via a discrete-time Riccati equation. Whenever possible, we consider the generalized discrete–time algebraic Riccati equation

0 = R(X) = C^{T} QC + A^{T} XA − E^{T} XE  (1)

−(A^{T} XB + C^{T} S)(R + B^{T} XB)^{−1}(B^{T} XA + S^{T} C),

where A, E ∈ ℝ^{n×n}, B ∈ ℝ^{n×m}, C ∈ ℝ^{p×n}, Q ∈ ℝ^{p×p}, R ∈ ℝ^{m×m}, and S ∈ ℝ^{p×m}. Furthermore, Q and R are assumed to be symmetric and A and E are large and spare. For the particular application above, we have

A = [ \begin{matrix}0 & 1 \\ & ... & ... \\ & & 0 & 1 \\ & & & 0 \end{matrix} ] .

The function R(X) is a rational matrix function, R(X) = 0 defines a system of nonlinear equations. Newton’s method for the numerical solution of DAREs can be formulated as follows

for k = 0, 1, 2, …

1. K_{k} ← K(X_{k}) = (R + B^{T} X_{k}B)^{−1}(B^{T} X_{k}A + S^{T} C).

2. A_{k} ← A − BK_{k}.

3. R_{k} ← R(X_{k}).

4. Solve for N_{k} in the Stein equation

A_{k}^{T} N_{k}A_{k} − E^{T} N_{k}E = −R_{k}.  (2)

5. X_{k+1} ← X_{k} + N_{k}. end for The computational cost for this algorithm mainly depends upon the cost for the numerical solution of the Stein equation (??). This can be done using the Bartels–Stewart algorithm [?] or an extension to the case E ≠ I [?, ?, ?]. The Bartels-Stewart algorithm is the standard direct method for the solution of Stein equations of small to moderate size. This method requires the computation of a Schur decomposition, and thus is not appropriate for large scale problems. The cost for the solution of the Stein equation is ≈ 73n^{3} flops. Iterative schemes have been developed including the Smith method [?], the sign-function method [?], and the alternating direction implicit (ADI) iteration method [?]. Unfortunately, all of these methods compute the solution in dense form and hence require O(n^{2}) storage. In case the solution to the Stein equation has low numerical rank (i.e., the eigenvalues decay rapidly) one can take advantage of this low rank structure to obtain approximate solutions in low rank factored form. If the effective rank is r ≪ n, then the storage is reduced from O(n^{2}) to O(nr). This approach will be discussed here in detail.

References

[1] R.H. Bartels and G.W. Stewart, Solution of the matrix equation AX + XB = C: Algorithm 432, Comm. ACM, 15 (1972), pp. 820–826.

[2] J.D. Gardiner, A.J. Laub, J.J. Amato, and C.B. Moler, Solution of the Sylvester matrix equation AXB + CXD = E, ACM Trans. Math. Software, 18 (1992), pp. 223–231.

[3] J.D. Gardiner, M.R. Wette, A.J. Laub, J.J. Amato, and C.B. Moler, Algorithm 705: A Fortran-77 software package for solving the Sylvester matrix equation AXB^{T} + CXD^{T} = E, ACM Trans. Math. Software, 18 (1992), pp. 232–238.

[4] T. Penzl, Numerical solution of generalized Lyapunov equations, Adv. Comp. Math., 8 (1997), pp. 33–48.

[5] J.D. Roberts, Linear model reduction and solution of the algebraic Riccati equation by use of the sign function, Internat. J. Control, 32 (1980), pp. 677–687. (Reprint of Technical Report No. TR-13, CUED/B-Control, Cambridge University, Engineering Department, 1971).

[6] A.H. Sayed and T. Kailath, A survey of spectral factorization methods, Num. Lin. Alg. Appl., 8 (2001), pp. 467–496.

[7] R.A. Smith, Matrix equation XA + BX = C, SIAM J. Appl. Math., 16 (1968), pp. 198–201.

[8] E.L. Wachspress, Iterative solution of the Lyapunov matrix equation, Appl. Math. Letters, 107 (1988), pp. 87–90.