Fast positive-real balanced truncation via quadratic alternating direction implicit iteration

Balanced truncation (BT), as applied to date in model order reduction (MOR), is known for its superior accuracy and computable error bounds. Positive-real BT (PRBT) is a particular BT procedure that preserves passivity and stability and imposes no structural constraints on the original state space. However, PRBT requires solving two algebraic Riccati equations (AREs), whose computational complexity limits its practical use in large-scale systems. This paper introduces a novel quadratic extension of the alternating direction implicit (ADI) iteration, which is called quadratic ADI (QADI), that efficiently solves an ARE. A Cholesky factor version of QADI, which is called CEQADI, exploits low-rank matrices and further accelerates PRBT. © 2007 IEEE.published_or_final_versio

Iterative and doubling algorithms for Riccati-type matrix equations: a comparative introduction

We review a family of algorithms for Lyapunov- and Riccati-type equations which are all related to each other by the idea of \emph{doubling}: they construct the iterate $Q_k = X_{2^k}$ of another naturally-arising fixed-point iteration $(X_h)$ via a sort of repeated squaring. The equations we consider are Stein equations $X - A^*XA=Q$, Lyapunov equations $A^*X+XA+Q=0$, discrete-time algebraic Riccati equations $X=Q+A^*X(I+GX)^{-1}A$, continuous-time algebraic Riccati equations $Q+A^*X+XA-XGX=0$, palindromic quadratic matrix equations $A+QY+A^*Y^2=0$, and nonlinear matrix equations $X+A^*X^{-1}A=Q$. We draw comparisons among these algorithms, highlight the connections between them and to other algorithms such as subspace iteration, and discuss open issues in their theory.Comment: Review article for GAMM Mitteilunge

Factorized solution of generalized stable Sylvester equations using many-core GPU accelerators

[EN] We investigate the factorized solution of generalized stable Sylvester equations such as those arising in model reduction, image restoration, and observer design. Our algorithms, based on the matrix sign function, take advantage of the current trend to integrate high performance graphics accelerators (also known as GPUs) in computer systems. As a result, our realisations provide a valuable tool to solve large-scale problems on a variety of platforms.We acknowledge support of the ANII - MPG Independent Research Group: "Efficient Hetergenous Computing" at UdelaR, a partner group of the Max Planck Institute in Magdeburg.Benner, P.; Dufrechou, E.; Ezzatti, P.; Gallardo, R.; Quintana-Ortí, ES. (2021). Factorized solution of generalized stable Sylvester equations using many-core GPU accelerators. The Journal of Supercomputing (Online). 77(9):10152-19164. https://doi.org/10.1007/s11227-021-03658-y101521916477

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

An inverse-free ADI algorithm for computing Lagrangian invariant subspaces

Summary: The numerical computation of Lagrangian invariant subspaces of large-scale Hamiltonian matrices is discussed in the context of the solution of Lyapunov equations. A new version of the low-rank alternating direction implicit method is introduced, which, in order to avoid numerical difficulties with solutions that are of very large norm, uses an inverse-free representation of the subspace and avoids inverses of ill-conditioned matrices. It is shown that this prevents large growth of the elements of the solution that may destroy a low-rank approximation of the solution. A partial error analysis is presented, and the behavior of the method is demonstrated via several numerical examples. Copyrigh

An inverse-free ADI algorithm for computing Lagrangian invariant subspaces

The numerical computation of Lagrangian invariant subspaces of large scale Hamiltonian matrices is discussed in the context of the solution of Lyapunov and Riccati equations. A new version of the low-rank alternating direction implicit method is introduced, which in order to avoid numerical difficulties with solutions that are of very large norm, uses an inverse-free representation of the subspace and avoids inverses of ill-conditioned matrices. It is shown that this prevents large growth of the elements of the solution which may destroy a low-rank approximation of the solution. A partial error analysis is presented and the behavior of the method is demonstrated via several numerical examples