A collection of functions has been implemented as a MatrixEquation.jl Julia package to solve several classes of Lyapunov, Sylvester and Riccati matrix equations. The goal was to demonstrate that programs written in the Julia language can achieve high computational performance, which is comparable with the performance of efficient structure exploiting Fortran implementations, as those available in the Systems and Control Library SLICOT.
The available functions cover both standard and generalized continuous and discrete Lyapunov, Sylvester and Riccati equations for both real and complex data. The functions for the solution of Lyapunov and Sylvester equations rely on efficient structure exploiting solvers for which the input data are in Schur or generalized Schur forms. A comprehensive set of Lyapunov and Sylvester operators has been implemented, which allow the estimation of condition numbers of these operators. The implementations of Riccati equation solvers employ orthogonal Schur vectors based methods and their extensions to linear matrix pencil based reduction approaches. The calls of all functions with adjoint (in complex case) or transposed (in real case) arguments are fully supported by appropriate computational algorithms, thus the matrix copying operations are mostly avoided.
The current version of the package includes the following main functions:
Solution of Lyapunov equations
lyapc Solution of the continuous Lyapunov equations AX+XA'+C = 0 and A'X+XA+C = 0 and of their generalized versions AXE'+EXA'+C = 0 and A'XE+E'XA+C = 0.
lyapd Solution of discrete Lyapunov equations AXA' - X +C = 0 and AXA' - X +C = 0 and of their generalized versions AXA' - EXE'+C = 0 and A'XA - E'XE+C = 0.
plyapc Computation of the upper triangular factor U of the positive semi-definite solution X = UU' of the continuous Lyapunov equation AX+XA'+BB' = 0 or of the positive semi-definite solution X = U'U of the continuous Lyapunov equation A'X+XA+B'B = 0 and of their generalized versions AXE' + EXA' + BB' = 0 and A'XE + E'XA + B'B = 0.
plyapd Computation of the upper triangular factor U of the positive semi-definite solution X = UU' of the discrete Lyapunov equation AXA' - X+BB' = 0 or of the positive semi-definite solution X = U'U of the discrete Lyapunov equation A'XA - X+B'B = 0 and of their generalized versions AXA' - EXE' + BB' = 0 and A'XA - E'XE + B'B = 0.
plyaps Efficient computation of the above solutions for A in a Schur form or the pair (A,E) in a generalized Schur form.
Solution of algebraic Riccati equations
arec Solution of the continuous Riccati equation A'X+XA-XRX+Q = 0 and of the control oriented continuous Riccati equation A'X+XA-(XB+S)R-1(B'X+S')+Q = 0.
garec Solution of the generalized continuous Riccati equation A'XE+E'XA-(A'XB+S)R-1(B'XA+S')+Q = 0.
ared Solution of the discrete Riccati equation A'XA - X - (A'XB+S)(R+B'XB)-1(B'XA+S') + Q = 0.
gared Solution of the generalized discrete Riccati equation A'XA - E'XE - (A'XB+S)(R+B'XB)-1(B'XA+S') + Q = 0.
Solution of Sylvester equations and systems
sylvc Solution of the (continuous) Sylvester equation AX+XB = C.
sylvd Solution of the (discrete) Sylvester equation AXB+X = C.
gsylv Solution of the generalized Sylvester equation AXB+CXD = E.
sylvsys Solution of the Sylvester system of matrix equations AX+YB = C, DX+YE = F.
dsylvsys Solution of the dual Sylvester system of matrix equations AX+DY = C, XB+YE = F.
Norm, condition and separation estimation of linear operators
opnorm1 Computation of the 1-norm of a linear operator.
opnorm1est Estimation of 1-norm of a linear operator.
oprcondest Estimation of the reciprocal 1-norm condition number of a linear operator.
opsepest Estimation of the separation of a linear operator.
A comprehensive set of Lyapunov and Sylvester operators and their inverses has been implemented. Each operator M allows to efficiently compute matrix-vector mutiplications M*x and M'*x. The inverse operators rely on the implemented Lyapunov and Sylvester equation solvers. All inverse operators have been implemented for both general matrices and for matrices in Schur or generalized Schur forms. This allows the efficient estimation of condition numbers and separations. The main solvers of Lyapunov and Sylvester equations rely on a set of specialized solvers for real or complex matrices in appropriate Schur forms. In what follows, for a matrix M, op(M) stays either for op(M) = M or op(M) = M'.
Solution of Lyapunov equations
lyapcs! Solution of the reduced continuous Lyapunov equations op(A)X + Xop(A)' + C = 0, with A in a Schur form, and op(A)Xop(E)' + op(E)Xop(A)' + C = 0, with (A,E) in a generalized Schur form.
lyapds! Solution of the reduced discrete Lyapunov equations op(A)Xop(A)' -X + C = 0, with A in a Schur form, and op(A)Xop(A)' - op(E)Xop(E)' + C = 0, with (A,E) in a generalized Schur form.
plyapcs! Computation of the upper triangular factor U of the positive solution X = op(U)op(U)' of the reduced continuous Lyapunov equation op(A)X + Xop(A)' + op(B)op(B)' = 0, with A in a Schur form, and of its generalized version op(A)Xop(E)' + op(E)Xop(A)' + op(B)op(B)' = 0, with (A,E) in a generalized Schur form.
plyapds! Computation of the upper triangular factor U of the positive solution X = op(U)op(U)' of the reduced discrete Lyapunov equation op(A)X*op(A)' - X + op(B)op(B)'=0, with A in a Schur form, and of its generalized version op(A)Xop(A)' + op(E)Xop(E)' + op(B)op(B)' = 0, with (A,E) in a generalized Schur form.
Solution of Sylvester equations and systems
sylvcs! Solution of the (continuous) Sylvester equation op1(A)X + Xop2(B) = C, with A and B in Schur form.
sylvds! Solution of the (discrete) Sylvester equation op1(A)Xop2(B) + X = C, with A and B in Schur form.
gsylvs! Solution of the generalized Sylvester equation op1(A)Xop2(B) + op1(C)Xop2(D) = E, with (A,C) and (B,D) in generalized Schur form.
sylvsyss! Solution of the Sylvester system of matrix equations AX+YB = C, DX+YE = F, with (A,C) and (B,D) in generalized Schur form.
dsylvsyss! Solution of the dual Sylvester system of matrix equations A'X+D'Y = C, XB'+YE' = F, with (A,C) and (B,D) in generalized Schur form.
tgsyl! LAPACK wrapper to solve the Sylvester systems of matrix equations AX-YB = C, DX-YE = F and A'X+D'Y = C, XB'+YE' = -F, with (A,D) and (B,E) in generalized Schur form.
For testing purposes, a set of solvers for Sylvester equations has been implemented, which employ the Kronecker-product expansion of the equations. These solvers are not recommended for large order matrices.
sylvckr Solution of the (continuous) Sylvester equation AX+XB = C.
sylvdkr Solution of the (discrete) Sylvester equation AXB+X = C.
gsylvkr Solution of the generalized Sylvester equation AXB+CXD = E.
sylvsyskr Solution of the Sylvester system of matrix equations AX+YB = C, DX+YE = F.
dsylvsyskr Solution of the dual Sylvester system of matrix equations AX+DY = C, XB+YE = F.