On the numerical solution of a class of systems of linear matrix equations

We consider the solution of systems of linear matrix equations in two or three unknown matrices. For dense problems we derive algorithms that determine the numerical solution by only involving matrices of the same size as those in the original problem, thus requiring low computational resources. For large and structured systems we show how the problem properties can be exploited to design effective algorithms with low memory and operation requirements. Numerical experiments illustrate the performance of the new methods

New optimization methods in predictive control

This thesis is mainly concerned with the efficient solution of a linear discrete-time finite horizon optimal control problem (FHOCP) with quadratic cost and linear constraints on the states and inputs. In predictive control, such a FHOCP needs to be solved online at each sampling instant. In order to solve such a FHOCP, it is necessary to solve a quadratic programming (QP) problem. Interior point methods (IPMs) have proven to be an efficient way of solving quadratic programming problems. A linear system of equations needs to be solved in each iteration of an IPM. The ill-conditioning of this linear system in the later iterations of the IPM prevents the use of an iterative method in solving the linear system due to a very slow rate of convergence; in some cases the solution never reaches the desired accuracy. A new well-conditioned IPM, which increases the rate of convergence of the iterative method is proposed. The computational advantage is obtained by the use of an inexact Newton method along with the use of novel preconditioners. A new warm-start strategy is also presented to solve a QP with an interior-point method whose data is slightly perturbed from the previous QP. The effectiveness of this warm-start strategy is demonstrated on a number of available online benchmark problems. Numerical results indicate that the proposed technique depends upon the size of perturbation and it leads to a reduction of 30-74% in floating point operations compared to a cold-start interior point method. Following the main theme of this thesis, which is to improve the computational efficiency of an algorithm, an efficient algorithm for solving the coupled Sylvester equation that arises in converting a system of linear differential-algebraic equations (DAEs) to ordinary differential equations is also presented. A significant computational advantage is obtained by exploiting the structure of the involved matrices. The proposed algorithm removes the need to solve a standard Sylvester equation or to invert a matrix. The improved performance of this new method over existing techniques is demonstrated by comparing the number of floating-point operations and via numerical examples

Least Squares Based Iterative Algorithm for the Coupled Sylvester Matrix Equations

By analyzing the eigenvalues of the related matrices, the convergence analysis of the least squares based iteration is given for solving the coupled Sylvester equations AX+YB=C and DX+YE=F in this paper. The analysis shows that the optimal convergence factor of this iterative algorithm is 1. In addition, the proposed iterative algorithm can solve the generalized Sylvester equation AXB+CXD=F. The analysis demonstrates that if the matrix equation has a unique solution then the least squares based iterative solution converges to the exact solution for any initial values. A numerical example illustrates the effectiveness of the proposed algorithm

Robust isogeometric preconditioners for the Stokes system based on the Fast Diagonalization method

In this paper we propose a new class of preconditioners for the isogeometric discretization of the Stokes system. Their application involves the solution of a Sylvester-like equation, which can be done efficiently thanks to the Fast Diagonalization method. These preconditioners are robust with respect to both the spline degree and mesh size. By incorporating information on the geometry parametrization and equation coefficients, we maintain efficiency on non-trivial computational domains and for variable kinematic viscosity. In our numerical tests we compare to a standard approach, showing that the overall iterative solver based on our preconditioners is significantly faster.Comment: 31 pages, 4 figure

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

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

Model reduction of controlled Fokker--Planck and Liouville-von Neumann equations

Model reduction methods for bilinear control systems are compared by means of practical examples of Liouville-von Neumann and Fokker--Planck type. Methods based on balancing generalized system Gramians and on minimizing an H2-type cost functional are considered. The focus is on the numerical implementation and a thorough comparison of the methods. Structure and stability preservation are investigated, and the competitiveness of the approaches is shown for practically relevant, large-scale examples

Coherent Quantum Filtering for Physically Realizable Linear Quantum Plants

The paper is concerned with a problem of coherent (measurement-free) filtering for physically realizable (PR) linear quantum plants. The state variables of such systems satisfy canonical commutation relations and are governed by linear quantum stochastic differential equations, dynamically equivalent to those of an open quantum harmonic oscillator. The problem is to design another PR quantum system, connected unilaterally to the output of the plant and playing the role of a quantum filter, so as to minimize a mean square discrepancy between the dynamic variables of the plant and the output of the filter. This coherent quantum filtering (CQF) formulation is a simplified feedback-free version of the coherent quantum LQG control problem which remains open despite recent studies. The CQF problem is transformed into a constrained covariance control problem which is treated by using the Frechet differentiation of an appropriate Lagrange function with respect to the matrices of the filter.Comment: 14 pages, 1 figure, submitted to ECC 201