A numerical comparison of solvers for large-scale, continuous-time algebraic Riccati equations and LQR problems

In this paper, we discuss numerical methods for solving large-scale continuous-time algebraic Riccati equations. These methods have been the focus of intensive research in recent years, and significant progress has been made in both the theoretical understanding and efficient implementation of various competing algorithms. There are several goals of this manuscript: first, to gather in one place an overview of different approaches for solving large-scale Riccati equations, and to point to the recent advances in each of them. Second, to analyze and compare the main computational ingredients of these algorithms, to detect their strong points and their potential bottlenecks. And finally, to compare the effective implementations of all methods on a set of relevant benchmark examples, giving an indication of their relative performance

On an integrated Krylov-ADI solver for large-scale Lyapunov equations

One of the most computationally expensive steps of the low-rank ADI method for large-scale Lyapunov equations is the solution of a shifted linear system at each iteration. We propose the use of the extended Krylov subspace method for this task. In particular, we illustrate how a single approximation space can be constructed to solve all the shifted linear systems needed to achieve a prescribed accuracy in terms of Lyapunov residual norm. Moreover, we show how to fully merge the two iterative procedures in order to obtain a novel, efcient implementation of the low-rank ADI method, for an important class of equations. Many state-of-the-art algorithms for the shift computation can be easily incorporated into our new scheme, as well. Several numerical results illustrate the potential of our novel procedure when compared to an implementation of the low-rank ADI method based on sparse direct solvers for the shifted linear systems

H2\mathcal{H}_2 Pseudo-Optimal Reduction of Structured DAEs by Rational Interpolation

In this contribution, we extend the concept of $\mathcal{H}_2$ inner product and $\mathcal{H}_2$ pseudo-optimality to dynamical systems modeled by differential-algebraic equations (DAEs). To this end, we derive projected Sylvester equations that characterize the $\mathcal{H}_2$ inner product in terms of the matrices of the DAE realization. Using this result, we extend the $\mathcal{H}_2$ pseudo-optimal rational Krylov algorithm for ordinary differential equations to the DAE case. This algorithm computes the globally optimal reduced-order model for a given subspace of $\mathcal{H}_2$ defined by poles and input residual directions. Necessary and sufficient conditions for $\mathcal{H}_2$ pseudo-optimality are derived using the new formulation of the $\mathcal{H}_2$ inner product in terms of tangential interpolation conditions. Based on these conditions, the cumulative reduction procedure combined with the adaptive rational Krylov algorithm, known as CUREd SPARK, is extended to DAEs. Important properties of this procedure are that it guarantees stability preservation and adaptively selects interpolation frequencies and reduced order. Numerical examples are used to illustrate the theoretical discussion. Even though the results apply in theory to general DAEs, special structures will be exploited for numerically efficient implementations

On a family of low-rank algorithms for large-scale algebraic Riccati equations

In [3] it was shown that four seemingly different algorithms for computing low-rank approximate solutions $X_j$ to the solution $X$ of large-scale continuous-time algebraic Riccati equations (CAREs) $0 = \mathcal{R}(X) := A^HX+XA+C^HC-XBB^HX$ generate the same sequence $X_j$ when used with the same parameters. The Hermitian low-rank approximations $X_j$ are of the form $X_j = Z_jY_jZ_j^H,$ where $Z_j$ is a matrix with only few columns and $Y_j$ is a small square Hermitian matrix. Each $X_j$ generates a low-rank Riccati residual $\mathcal{R}(X_j)$ such that the norm of the residual can be evaluated easily allowing for an efficient termination criterion. Here a new family of methods to generate such low-rank approximate solutions $X_j$ of CAREs is proposed. Each member of this family of algorithms proposed here generates the same sequence of $X_j$ as the four previously known algorithms. The approach is based on a block rational Arnoldi decomposition and an associated block rational Krylov subspace spanned by $A^H$ and $C^H.$ Two specific versions of the general algorithm will be considered; one will turn out to be a rediscovery of the RADI algorithm, the other one allows for a slightly more efficient implementation compared to the RADI algorithm. Moreover, our approach allows for adding more than one shift at a time

Effizientes Lösen von großskaligen Riccati-Gleichungen und ein ODE-Framework für lineare Matrixgleichungen

This work considers the iterative solution of large-scale matrix equations. Due to the size of the system matrices in large-scale Riccati equations the solution can not be calculated directly but is approximated by a low rank matrix ZYZ^*. Herein Z is a basis of a low-dimensional rational Krylov subspace. The inner matrix Y is a small square matrix. Two ways to choose this inner matrix are examined: By imposing a rank condition on the Riccati residual and by projecting the Riccati residual onto the Krylov subspace generated by Z. The rank condition is motivated by the well-known ADI iteration. The ADI solutions span a rational Krylov subspace and yield a rank-p residual. It is proven that the rank-p condition guarantees existence and uniqueness of such an approximate solution. Known projection methods are generalized to oblique projections and a new formulation of the Riccati residual is derived, which allows for an efficient evaluation of the residual norm. Further a truncated approximate solution is characterized as the solution of a Riccati equation, which is projected to a subspace of the Krylov subspace generated by Z. For the approximate solution of Lyapunov equations a system of ordinary differential equations (ODEs) is solved via Runge-Kutta methods. It is shown that the space spanned by the approximate solution is a rational Krylov subspace with poles determined by the time step sizes and the eigenvalues of the matrices of the Butcher tableau of the used Runge-Kutta method. The method is applied to a model order reduction problem. The analytical solution of the system of ODEs satisfies an algebraic invariant. Those Runge-Kutta methods which preserve this algebraic invariant are characterized by a simple condition on the corresponding Butcher tableau. It is proven that these methods are equivalent to the ADI iteration. The invariance approach is transferred to Sylvester equations.Diese Arbeit befasst sich mit der numerischen Lösung hochdimensionaler Matrixgleichungen mittels iterativer Verfahren. Aufgrund der Größe der Systemmatrizen in großskaligen algebraischen Riccati-Gleichung kann die Lösung nicht direkt bestimmt werden, sondern wird durch eine approximative Lösung ZYZ^* von geringem Rang angenähert. Hierbei wird Z als Basis eines rationalen Krylovraums gewählt und enthält nur wenige Spalten. Die innere Matrix Y ist klein und quadratisch. Es werden zwei Wege untersucht, die Matrix Y zu wählen: Durch eine Rang-Bedingung an das Riccati-Residuum und durch Projektion des Riccati-Residuums auf den von Z erzeugten Krylovraum. Die Rang-Bedingung wird durch die wohlbekannten ADI-Verfahren motiviert. Die approximativen ADI-Lösungen spannen einen Krylovraum auf und führen zu einem Riccati-Residuum vom Rang p. Es wird bewiesen, dass die Rang-p-Bedingung Existenz und Eindeutigkeit einer solchen approximativen Lösung impliziert. Aus diesem Ergebnis werden effiziente iterative Verfahren abgeleitet, die eine solche approximative Lösung erzeugen. Bisher bekannte Projektionsverfahren werden auf schiefe Projektionen erweitert und es wird eine neue Formulierung des Riccati-Residuums hergeleitet, die eine effiziente Berechnung der Norm erlaubt. Weiter wird eine abgeschnittene approximative Lösung als Lösung einer Riccati-Gleichung charakterisiert, die auf einen Unterraum des von Z erzeugten Krylovraums projiziert wird. Um die Lösung der Lyapunov-Gleichung zu approximieren wird ein System gewöhnlicher Differentialgleichungen mittels Runge-Kutta-Verfahren numerisch gelöst. Es wird gezeigt, dass der von der approximativen Lösung aufgespannte Raum ein rationaler Krylovraum ist, dessen Pole von den Zeitschrittweiten der Integration und den Eigenwerten der Koeffizientenmatrix aus dem Butcher-Tableau des verwendeten Runge-Kutta-Verfahrens abhängen. Das Verfahren wird auf ein Problem der Modellreduktion angewendet. Die analytische Lösung des Differentialgleichungssystems erfüllt eine algebraische Invariante. Diejenigen Runge-Kutta-Verfahren, die diese Invariante erhalten, werden durch eine Bedingung an die zugehörigen Butcher-Tableaus charakterisiert. Es wird gezeigt, dass diese speziellen Verfahren äquivalent zur ADI-Iteration sind. Der Invarianten-Ansatz wird auf Sylvester-Gleichungen übertragen