1. Federico Poloni
  2. pgdoubling




a MATLAB package to solve algebraic Riccati equations and optimal control problems using permuted graph bases

© Federico Poloni, 2012 -- 2014


This package uses the algorithms described in [MehP12, MehP13] to compute Lagrangian invariant subspaces of Hamiltonian and "extended Hamiltonian" pencils. The resulting algorithm is more accurate than Matlab's care on several tricky test problems [carex].

It is the right choice for you if:

  • you want to compute a linear-quadratic regulator, solve a Riccati equation, or a generalized Riccati equation (those that correspond to a matrix pencil in the form:
    [ 0  I 0]    [0  A  B]
    [ -I 0 0]s - [A' Q  S]
    [ 0  0 0]    [B' S' R]


  • you are not satisfied with the accuracy obtained by Matlab's care, or wish to test a new algorithm;

  • your problem is small-scale (up to a few hundred degrees of freedom);

  • your system is not a descriptor system (i.e., no E coefficient, just \dot{x} = Ax + Bu).

  • Note: your problem could also be a singular control problem (Lur'e equations, singular R); there is also a function to compute the optimal parameter gamma in H_infinity control;

The current implementation is in pure Matlab (no mex-files), so it is slower that care. Rewriting the more intensive parts in Fortran or C will very likely result in a huge speedup.


The files named test*.m are tests. In order to run them, you need to install and have in your Matlab path (change it using pathtool) the following two libraries:

The command runtests should take about one minute and return success. You will see some warnings about ill-conditioned matrices being inverted, but it is ok.


A CanBasis is a basis matrix for an n-dimensional subspace in the form P*[I(n);X], where P is a permutation matrix and X is arbitrary. Usually the routines return CanBases in which each element of X is smaller (in modulus) than a certain configurable threshold. A CanBasis is stored as a struct can with fields can.X and can.p (permutation vector).

A SymBasis is the structured version of a CanBasis for a Lagrangian subspace. P is a "symplectic swap matrix" (a symplectic analogous of a permutation), and X is Hermitian. A SymBasis is stored as a struct sym with fields sym.X and sym.v (boolean vector).

You can represent subspaces, matrices and matrix pencils through canBases and symBases. See the basic theory in [MehP12]. (all will be clearer if you have read this paper.)


All functions contain help; type help functionName in Matlab to access it. You can look in the "test" functions for some more usage examples. In the following there is a quick introduction.

Basic usage: matrix equations


Returns the symmetric (semi)stabilizing solution of a CARE A'X+XA+Q=XGX. Also returns the (semi)stabilizing solution of the "dual CARE" YA'+AY+YQY=G. More important, U and V are basis matrices for the (semi)stabilizing and (semi)antistabilizing Lagrangian subspaces of the Hamiltonian, i.e., subspaces U=[U1;U2] and V=[V1;V2] such that X=U2/U1, Y=V1/V2. You should always try to restate your successive computations in order to use U and V rather than X and Y, as X and Y may be ill-conditioned while U and V are tame.

Recommended options for robustness: solveCARE(A,G,Q,'type','sign','safer',true). If your problem has large size, you may consider increasing the threshold with (...'threshold',sqrt(n),'diagonalThreshold',sqrt(n)). The impact of the thresholds hasn't been investigated thoroughly yet.


solve a (possibly singular) generic control problem in the form:

    [ 0  I 0]    [0  A  B]
    [ -I 0 0]s - [A' Q  S],
    [ 0  0 0]    [B' S' R]

returning the (semi)stabilizing X and Y. In addtion, returns bases U and V for the stabilizing subspace [X;I] and the antistabilizing one [I;Y]. Note that the position of X and I in the subspace is different from the one in the 2x2 problem. I am just following the standard notation for each of the two problems. In any case, you should work with U and V rather than with X and Y, as argued above.


run the gamma iteration [BenBMX07,MehP13] on a H_inf control problem, with tolerance tol.

Some more matrix equations could be added --- These techniques work also for DARE, X=Q+A*inv(X)A', X=Q-Ainv(X)*A', and 0=P+QY+RY\^2 (with some modifications). If you would like to solve one of these with this package, contact the author.

Basic usage: building blocks

canBasisFromSubspace takes a subspace and returns X and a permutation P defining its CanBasis. symBasisFromSymplecticSubspace takes a Lagrangian subspace and returns its SymBasis. The associated permutation matrix P is obtained by rowSwap(eye(2*n),sym.v,'N'), if needed.

You may also be interested in symBasisFromSymplecticPencil, symBasisFromHamiltonianPencil, symBasisFromEvenPencil (actually only works for even pencils in the form above; a generic even pencil does not have a symBasis).

Some random stuff that an end user might find interesting: checkCAREInvariantSubspaceResidual (computes and returns several error measures), riccatiCoefficientsFromHamiltonian, hamiltonian (should be named hamiltonianFromRiccatiCoefficients), pirq ((symplectic swap P)RQ factorization), ChuLM07Carex (returns the results of the experiments shown in [ChuLM07]), randomLagrangianSubspace.


This library does not have a large user base; I am happy to hear about end users' experiences ("I tried it and it works great; I tried it and it doesn't work") and to receive feedback and bug reports.

If you use it in a paper, please cite us when it is relevant. In the academia we basically live off citation counts.


Check COPYING.txt to see if you qualify to use, modify and redistribute the library in other projects for free.


Q: why are the function names so long?

A: I find it more expressive than the default 'shorten everything' Matlab syntax. You can configure Matlab's editor for tab completion and the fact that the names are long will not slow you down.

Q: why is it so slow?

A: mainly, because the PiRQ factorization is implemented with a loop in pure Matlab. Rewriting it in a MEX-file would result in a massive speed improvement on large problems. In general, the functions are not aggressively optimized; some more performance could be squeezed out.

Q: why do you roll your own option parsing functions? OOP is slow in Matlab!

A: I haven't found yet an option parser for Matlab that satisfies some basic requirements, namely: (1) can pass options unchanged to subroutines (2) complains at the end if there are unused parameters (3) has an acceptable syntax which does not require me to specify option names in 18 different places. If you know of one, let me know.

Q: does it work when there is a mass matrix E?

A: no, and it would require significant changes to the algorithm to support it. We'll be working on that hopefully.


Mehrmann, Volker; Poloni, Federico Doubling algorithms with permuted Lagrangian graph bases. SIAM J. Matrix Anal. Appl. 33 (2012), no. 3, 780–805.
Mehrmann, Volker; Poloni, Federico Using permuted graph bases in H∞ control. Automatica J. IFAC 49 (2013), no. 6, 1790–1797.
Chu, Delin; Liu, Xinmin; Mehrmann, Volker A numerical method for computing the Hamiltonian Schur form. Numer. Math. 105 (2007), no. 3, 375–412.
Benner, Peter; Byers, Ralph; Mehrmann, Volker; Xu, Hongguo A robust numerical method for the γ-iteration in H∞ control. Linear Algebra Appl. 425 (2007), no. 2-3, 548–570.
Jörn Abels , Peter Benner CAREX - A Collection of Benchmark Examples for Continuous-Time Algebraic Riccati Equations (Version 2.0) (1999), http://citeseerx.ist.psu.edu/viewdoc/summary?doi=