Riccati and Quadratic Matrix Equation Solvers

arec(A, G, Q = 0; scaling = 'B', pow2 = false, as = false, rtol::Real = nϵ, nrm = 1) -> (X, EVALS, Z, scalinfo)

Compute X, the hermitian/symmetric stabilizing solution (if as = false) or anti-stabilizing solution (if as = true) of the continuous-time algebraic Riccati equation

 A'X + XA - XGX + Q = 0,

where G and Q are hermitian/symmetric matrices or uniform scaling operators. Scalar-valued G and Q are interpreted as appropriately sized uniform scaling operators G*I and Q*I. The Schur method of [1] is used.

To enhance the accuracy of computations, a block scaling of matrices G and Q is performed, if the default setting scaling = 'B' is used. This scaling is however performed only if norm(Q) > norm(G). A general, eigenvalue computation oriented scaling combined with a block scaling is used if scaling = 'G' is selected. An alternative, experimental structure preserving scaling can be performed using the option scaling = 'S'. A symmetric matrix equilibration based scaling is employed if scaling = 'K', for which the underlying vector norm can be specified using the keyword argument nrm = p, where p = 1 is the default setting. Scaling can be disabled with the choice scaling = 'N'. If pow2 = true, the scaling elements are enforced to the nearest power of 2 (default: pow2 = false).

By default, the lower bound for the 1-norm reciprocal condition number rtol is n*ϵ, where n is the order of A and ϵ is the machine epsilon of the element type of A.

EVALS is a vector containing the (stable or anti-stable) eigenvalues of A-GX.

Z = [U; V] is an orthogonal basis for the stable/anti-stable deflating subspace such that X = Sx*(V/U)*Sxi, where Sx and Sxi are diagonal scaling matrices contained in the named tuple scalinfo as scalinfo.Sx and scalinfo.Sxi, respectively.

Note: To solve the continuous-time algebraic Riccati equation

 A'X + XA - XBR^(-1)B'X + Q = 0,

with R a hermitian/symmetric matrix and B a compatible size matrix, G = BR^(-1)B' must be provided. This approach is not numerically suited when R is ill-conditioned and/or B has large norm.

Reference:

[1] Laub, A.J., A Schur Method for Solving Algebraic Riccati equations. IEEE Trans. Auto. Contr., AC-24, pp. 913-921, 1979.

Example

julia> using LinearAlgebra

julia> A = [-6. -2. 1.; 5. 1. -1; -4. -2. -1.]
3×3 Array{Float64,2}:
 -6.0  -2.0   1.0
  5.0   1.0  -1.0
 -4.0  -2.0  -1.0

julia> G = [1. 0. 0.; 0. 5. 0.; 0. 0. 10.]
3×3 Array{Float64,2}:
 1.0  0.0   0.0
 0.0  5.0   0.0
 0.0  0.0  10.0

julia> X, CLSEIG = arec(A,G,2I);

julia> X
3×3 Array{Float64,2}:
  0.459589   0.333603   -0.144406
  0.333603   0.65916    -0.0999216
 -0.144406  -0.0999216   0.340483

julia> A'*X+X*A-X*G*X+2I
3×3 Array{Float64,2}:
  2.22045e-16  4.44089e-16  -1.77636e-15
  4.44089e-16  6.66134e-16   1.11022e-16
 -1.77636e-15  1.11022e-16  -1.33227e-15

julia> CLSEIG
3-element Array{Complex{Float64},1}:
 -4.411547592296008 + 2.4222082620381102im
 -4.411547592296008 - 2.4222082620381102im
 -4.337128244724371 + 0.0im

julia> eigvals(A-G*X)
3-element Array{Complex{Float64},1}:
 -4.4115475922960075 - 2.4222082620381076im
 -4.4115475922960075 + 2.4222082620381076im
  -4.337128244724374 + 0.0im

arec(A, B, R, Q, S; scaling = 'B', pow2 = false, as = false, rtol::Real = nϵ, orth = false, nrm = 1) -> (X, EVALS, F, Z, scalinfo)

Compute X, the hermitian/symmetric stabilizing solution (if as = false) or anti-stabilizing solution (if as = true) of the continuous-time algebraic Riccati equation

 A'X + XA - (XB+S)R^(-1)(B'X+S') + Q = 0,

where R and Q are hermitian/symmetric matrices or uniform scaling operators such that R is nonsingular. Scalar-valued R and Q are interpreted as appropriately sized uniform scaling operators R*I and Q*I. S, if not specified, is set to S = zeros(size(B)). The Schur method of [1] is used.

To enhance the accuracy of computations, a block scaling of matrices R, Q and S is performed, if the default setting scaling = 'B' is used. This scaling is however performed only if norm(Q) > norm(B)^2/norm(R). A general, eigenvalue computation oriented scaling combined with a block scaling is used if scaling = 'G' is selected. An alternative, structure preserving scaling can be performed using the option scaling = 'S'. A symmetric matrix equilibration based scaling is employed if scaling = 'K', for which the underlying vector norm can be specified using the keyword argument nrm = p, where p = 1 is the default setting. Experimental structure preserving scalings can be performed using the options scaling = 'D' or scaling = 'T'. Scaling can be disabled with the choice scaling = 'N'. If pow2 = true, the scaling elements are enforced to the nearest power of 2 (default: pow2 = false).

By default, the lower bound for the 1-norm reciprocal condition number rtol is n*ϵ, where n is the order of A and ϵ is the machine epsilon of the element type of A.

EVALS is a vector containing the (stable or anti-stable) eigenvalues of A-BF.

F is the stabilizing or anti-stabilizing gain matrix F = R^(-1)(B'X+S').

Z = [U; V; W] is a basis for the relevant stable/anti-stable deflating subspace such that X = Sx*(V/U)*Sxi and F = -Sr*(W/U)*Sxi, where Sx, Sxi and Sr are diagonal scaling matrices contained in the named tuple scalinfo as scalinfo.Sx, scalinfo.Sxi and scalinfo.Sr, respectively. An orthogonal basis Z can be determined, with an increased computational cost, by setting orth = true.

Reference:

[1] Laub, A.J., A Schur Method for Solving Algebraic Riccati equations. IEEE Trans. Auto. Contr., AC-24, pp. 913-921, 1979.

Example

julia> using LinearAlgebra

julia> A = [-6. -2. 1.; 5. 1. -1; -4. -2. -1.]
3×3 Array{Float64,2}:
 -6.0  -2.0   1.0
  5.0   1.0  -1.0
 -4.0  -2.0  -1.0

julia> B = [1. 2.; 2. 0.; 0. 1.]
3×2 Array{Float64,2}:
 1.0  2.0
 2.0  0.0
 0.0  1.0

julia> R = [1. 0.; 0. 5.]
2×2 Array{Float64,2}:
 1.0  0.0
 0.0  5.0

julia> X, CLSEIG, F = arec(A,B,R,2I);

julia> X
3×3 Array{Float64,2}:
  0.522588   0.303007  -0.327227
  0.303007   0.650895  -0.132608
 -0.327227  -0.132608   0.629825

julia> A'*X+X*A-X*B*inv(R)*B'*X+2I
3×3 Array{Float64,2}:
 -2.66454e-15  -1.55431e-15   8.88178e-16
 -1.55431e-15   2.22045e-15  -2.9976e-15
  9.99201e-16  -2.9976e-15    4.44089e-16

julia> CLSEIG
3-element Array{Complex{Float64},1}:
   -4.37703628399912 + 2.8107164873731247im
   -4.37703628399912 - 2.8107164873731247im
 -1.8663764577096091 + 0.0im

julia> eigvals(A-B*F)
3-element Array{Complex{Float64},1}:
  -4.377036283999118 - 2.8107164873731234im
  -4.377036283999118 + 2.8107164873731234im
 -1.8663764577096063 + 0.0im

arec(A, B, G, R, Q, S; scaling = 'B', pow2 = false, as = false, rtol::Real = nϵ, orth = false, nrm = 1) -> (X, EVALS, F, Z, scalinfo)

Computes X, the hermitian/symmetric stabilizing solution (if as = false) or anti-stabilizing solution (if as = true) of the continuous-time algebraic Riccati equation

 A'X + XA - XGX - (XB+S)R^(-1)(B'X+S') + Q = 0,

where G, R and Q are hermitian/symmetric matrices or uniform scaling operators such that R is nonsingular. Scalar-valued G, R and Q are interpreted as appropriately sized uniform scaling operators G*I, R*I and Q*I. S, if not specified, is set to S = zeros(size(B)). For well conditioned R, the Schur method of [1] is used. For ill-conditioned R or if orth = true, the generalized Schur method of [2] is used.

To enhance the accuracy of computations, a block oriented scaling of matrices G, R, Q and S is performed using the default setting scaling = 'B'. This scaling is performed only if norm(Q) > max(norm(G), norm(B)^2/norm(R)). A general, eigenvalue computation oriented scaling combined with a block scaling is used if scaling = 'G' is selected. An alternative, structure preserving scaling can be performed using the option scaling = 'S'. A symmetric matrix equilibration based scaling is employed if scaling = 'K', for which the underlying vector norm can be specified using the keyword argument nrm = p, where p = 1 is the default setting. If orth = true, two experimental scaling procedures can be activated using the options scaling = 'D' and scaling = 'T'. Scaling can be disabled with the choice scaling = 'N'.

By default, the lower bound for the 1-norm reciprocal condition number rtol is n*ϵ, where n is the order of A and ϵ is the machine epsilon of the element type of A.

EVALS is a vector containing the (stable or anti-stable) eigenvalues of A-BF-GX.

F is the stabilizing or anti-stabilizing gain matrix F = R^(-1)(B'X+S').

Z = [U; V; W] is a basis for the relevant stable/anti-stable deflating subspace such that X = Sx*(V/U)*Sxi and F = -Sr*(W/U)*Sxi, where Sx, Sxi and Sr are diagonal scaling matrices contained in the named tuple scalinfo as scalinfo.Sx, scalinfo.Sxi and scalinfo.Sr, respectively. An orthogonal basis Z can be determined, with an increased computational cost, by setting orth = true.

Reference:

[1] Laub, A.J., A Schur Method for Solving Algebraic Riccati equations. IEEE Trans. Auto. Contr., AC-24, pp. 913-921, 1979.

[2] W.F. Arnold, III and A.J. Laub, Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations, Proc. IEEE, 72:1746-1754, 1984.

ared(A, B, R, Q, S; scaling = 'B', pow2 = false, as = false, rtol::Real = nϵ, nrm = 1) -> (X, EVALS, F, Z, scalinfo)

Compute X, the hermitian/symmetric stabilizing solution (if as = false) or anti-stabilizing solution (if as = true) of the discrete-time algebraic Riccati equation

A'XA - X - (A'XB+S)(R+B'XB)^(-1)(B'XA+S') + Q = 0,

where R and Q are hermitian/symmetric matrices. Scalar-valued R and Q are interpreted as appropriately sized uniform scaling operators R*I and Q*I. S, if not specified, is set to S = zeros(size(B)).

To enhance the accuracy of computations, a block oriented scaling of matrices R, Q and S is performed using the default setting scaling = 'B'. This scaling is performed only if norm(Q) > norm(B)^2/norm(R). A general, eigenvalue computation oriented scaling combined with a block scaling is used if scaling = 'G' is selected. An alternative, structure preserving scaling can be performed using the option scaling = 'S'. A symmetric matrix equilibration based scaling is employed if scaling = 'K', for which the underlying vector norm can be specified using the keyword argument nrm = p, where p = 1 is the default setting. Experimental structure preserving scalings can be performed using the options scaling = 'D', scaling = 'R' and scaling = 'T'. Scaling can be disabled with the choice scaling = 'N'. If pow2 = true, the scaling elements are enforced to the nearest power of 2 (default: pow2 = false).

By default, the lower bound for the 1-norm reciprocal condition number rtol is n*ϵ, where n is the order of A and ϵ is the machine epsilon of the element type of A.

EVALS is a vector containing the (stable) eigenvalues of A-BF.

F is the stabilizing gain matrix F = (R+B'XB)^(-1)(B'XA+S').

Z = [U; V; W] is an orthogonal basis for the relevant stable/anti-stable deflating subspace such that X = Sx*(V/U)*Sxi and F = -Sr*(W/U)*Sxi, where Sx, Sxi and Sr are diagonal scaling matrices contained in the named tuple scalinfo as scalinfo.Sx, scalinfo.Sxi and scalinfo.Sr, respectively.

Reference:

[1] W.F. Arnold, III and A.J. Laub, Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations, Proc. IEEE, 72:1746-1754, 1984.

Example

julia> using LinearAlgebra

julia> A = [ 0. 1.; 0. 0. ]
2×2 Array{Float64,2}:
 0.0  1.0
 0.0  0.0

julia> B = [ 0.; sqrt(2.) ]
2-element Array{Float64,1}:
 0.0
 1.4142135623730951

julia> R = 1.
1.0

julia> Q = [ 1. -1.; -1. 1. ]
2×2 Array{Float64,2}:
  1.0  -1.0
 -1.0   1.0

julia> X, CLSEIG, F = ared(A,B,R,Q);

julia> X
2×2 Array{Float64,2}:
  1.0  -1.0
 -1.0   1.5

julia> A'*X*A-X-A'*X*B*inv(R+B'*X*B)*B'*X*A+Q
2×2 Array{Float64,2}:
  0.0          -3.33067e-16
 -3.33067e-16   8.88178e-16

julia> CLSEIG
2-element Array{Complex{Float64},1}:
 0.4999999999999998 - 0.0im
               -0.0 - 0.0im

julia> eigvals(A-B*F)
2-element Array{Float64,1}:
 -2.7755575615628914e-16
  0.5

nare(A1, A2, G, Q; scaling = 'B', pow2 = false, as = false, disc = false, stab = true, 
     rtol::Real = nϵ, nrm = 1, indeig) -> (X, EVALS, Z, scalinfo, ind, HEVALS)

Compute X, the stabilizing solution (if as = false) or anti-stabilizing solution (if as = true) of the nonsymmetric/nonhermitian algebraic Riccati equation (NARE)

 A1X + XA2 - XGX + Q = 0,

where A1, A2, G and Q are m×m, n×n, n×m and m×n matrices, respectively, and X is an m×n matrix. An extension of the Schur-form technique of [1] is used, where a suitable invariant subspace [U; V] corresponding to the selected n eigenvalues of a Hamiltonian-like matrix H are determined and X = V*inv(U) is the resulting solution. Solvents corrsponding to arbitrary collections of n eigenvalues of H can be selected by specifying their indices via indeig = ind, where ind must be an integer vector with at least n uniques indices being part of a permutation of the vector 1:n+m. For real data, the indices must additionally indicate a symmetric set of complex conjugate eigenvalues.

To enhance the accuracy of computations, a block scaling of matrices G and Q is performed, if the default setting scaling = 'B' is used. This scaling is however performed only if norm(Q) > norm(G) and norm(G) > 0. A general, eigenvalue computation oriented scaling combined with a block scaling is used if scaling = 'G' is selected. If n = m, an alternative, experimental structure preserving scaling can be performed using the option scaling = 'S'. For the same case, a symmetric matrix equilibration based scaling can be employed if scaling = 'K', for which the underlying vector norm can be specified using the keyword argument nrm = p, where p = 1 is the default setting. Scaling can be disabled with the choice scaling = 'N'. If pow2 = true, the scaling elements are enforced to the nearest power of 2 (default: pow2 = false).

By default, the lower bound for the 1-norm reciprocal condition number rtol is n*ϵ, where n is the order of A and ϵ is the machine epsilon of the element type of A1.

X is a stabilizing or anti-stabilizing solvent of the Riccati equation for the selected/chosen eigenvalues of the Hamiltonian-like matrix H := [A1 -G; -Q -A2].0 -B2' 0 ].X = nothingif no such a solution exists, in which caseZ,indandHEVALS` (see below) contain additional information to serve for possible alternative choices.

EVALS is a vector containing n stable or n anti-stable eigenvalues of H, if stab = true. These are also the eigenvalues of A2-GX. If stab = false, the eigenvalues must not be all stable or all anti-stable. The selection of the eigenvalues is performed as follows. If disc = false, then the n eigenvalues with the smallest real parts are selected if as = false, or with the largest real parts if as = true. If disc = true, then the n eigenvalues with the smallest moduli are selected if as = false, or with the largest moduli if as = true. Possible permutations of eigenvalues can occur if the n-th eigenvalue belongs to a complex conjugated pair.

Z = [U; V] is an orthogonal basis for the invariant subspace corresponding to EVALS, such that X = Sx*(V/U)*Sxi, where Sx and Sxi are diagonal scaling matrices contained in the named tuple scalinfo as scalinfo.Sx and scalinfo.Sxi, respectively.

ind contains the indices of the eigenvalues used to compute the invariant subspace underlying the computation of the resulting solvent X.

HEVALS is a vector containing n+m eigenvalues of the Hamiltonian-like matrix H. Only HEVALS can be computed with the call HEVALS = nare(A1, A2, G, Q; onlyeigvals = true, ...).

Note 1: To solve the continuous-time nonsymmetric/nonhermitian algebraic Riccati equation

 A1'X + XA2 - XB1R^(-1)B2'X + Q = 0,

with R a nonsingular matrix and B1 and B2 of compatible size matrices, G = B1R^(-1)B2' must be provided. This approach is not numerically suited when R is ill-conditioned and/or B1 and B2 have large norms.

Note 2: Generically, for complex data, there are up most C(n+m,n) complex solvents which correspond to all possible combinations of n eigenvalues of the n+m eigenvalues of the matrix H. For real data, there are usually less solvents than for complex data. All solvents can be determined by generating an exhaustive set of C(n+m,n) permutations and checking for each case if a solution exists (see some of the attached tests examples using try/catch constructs).

Reference:

[1] Laub, A.J., A Schur Method for Solving Algebraic Riccati equations. IEEE Trans. Auto. Contr., AC-24, pp. 913-921, 1979.

narec(A1, A2, B1, B2, R, Q, S1 = 0, S2 = 0; scaling = 'B', pow2 = false, as = false, rtol::Real = nϵ, orth = false, nrm = 1, indeig) -> (X, EVALS, F, Z, scalinfo, ind, HEVALS)

Computes X, the stabilizing solution (if as = false) or anti-stabilizing solution (if as = true) of the continuous-time nonsymmetric/nonhermitian algebraic Riccati equation

 A2'X + XA1 - (XB1+S1)R^(-1)(B2'X+S2') + Q = 0,

where A1 is an n×n matrix, A2 is an l×l matrix, and B1, B2, R, Q, S1, S2 are n×m, l×m, m×m, l×n, l×m, n×m matrices, respectively, with R nonsingular. Scalar-valued R is interpreted as an appropriately sized uniform scaling operator R*I. S1 and S2, if not specified, are set to zero matrices of appropriate sizes. For well conditioned R, the Schur method of [1] is used for a suitably formulated nonsymmetrix/nonhermitian algebraic Riccati equation. For ill-conditioned R or if orth = true, the generalized Schur method of [2] and [3] is used. Solvents corrsponding to arbitrary collections of n generalized finite eigenvalues of the pair (H,J) (see below) can be selected by specifying their indices via indeig = ind, where ind must be an integer vector with at least n uniques indices being part of a permutation of the vector 1:n+l. For real data, the indices must additionally indicate a symmetric set of complex conjugate eigenvalues.

To enhance the accuracy of computations, a block oriented scaling of matrices R, Q and S1 is performed using the default setting scaling = 'B'. This scaling is performed only if norm(Q) > norm(B1)^2/norm(R). A general, eigenvalue computation oriented scaling combined with a block scaling is used if scaling = 'G' is selected. If A1 and A2 have the same sizes, an alternative, structure preserving scaling can be performed using the option scaling = 'S'. For the same case, a symmetric matrix equilibration based scaling can be employed if scaling = 'K', for which the underlying vector norm can be specified using the keyword argument nrm = p, where p = 1 is the default setting. If orth = true, two experimental scaling procedures can be activated using the options scaling = 'D' or scaling = 'T' (the later only if A1 and A2 have the same sizes). Scaling can be disabled with the choice scaling = 'N'.

By default, the lower bound for the 1-norm reciprocal condition number rtol is n*ϵ, where n is the order of A1 and ϵ is the machine epsilon of the element type of A1.

X is a stabilizing or anti-stabilizing solvent of the Riccati equation for the selected/chosen finite generalized eigenvalues of the pair (H,J), where H = [ A1 0 B1; -Q -A2' -S1; S2' B2' R ] and J = [ I 0 0; 0 I 0; 0 0 0 ]. X = nothing if no such a solution exists, in which case Z, ind and HEVALS (see below) contain additional information to serve for possible alternative choices.

EVALS is a vector containing n stable or n anti-stable finite generalized eigenvalues of the pair (H,J), which underlies the computation of the solution. These are also the eigenvalues of A1-B1F.

F is the stabilizing or anti-stabilizing gain matrix F = R^(-1)(B2'X+S2').

Z = [U; V; W] is a basis for the relevant stable/anti-stable deflating subspace such that X = Sx*(V/U)*Sxi and F = -Sr*(W/U)*Sxi, where Sx, Sxi and Sr are diagonal scaling matrices contained in the named tuple scalinfo as scalinfo.Sx, scalinfo.Sxi and scalinfo.Sr, respectively. An orthogonal basis Z can be determined with an increased computational cost, by setting orth = true.

ind contains the indices of the eigenvalues used to compute the deflating subspace underlying the computation of the resulting solvent X.

HEVALS is a vector containing the n+l finite eigenvalues of the pair (H,J). Only HEVALS can be computed with the call HEVALS = narec(A1, A2, B1, B2, R, Q, S1 = 0, S2 = 0; onlyeigvals = true, ...).

Note: The considered nonsymmetric/nonhermitian Riccati equation is related to the computation of Nash equilibria in non-cooperative differential games [4]. The option to select eigenvalues via indeig = ind is instrumental for determining various solvents.

Reference:

[1] Laub, A.J., A Schur Method for Solving Algebraic Riccati equations. IEEE Trans. Auto. Contr., AC-24, pp. 913-921, 1979.

[2] W.F. Arnold, III and A.J. Laub, Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations, Proc. IEEE, 72:1746-1754, 1984.

[3] M. Jungers and C. Oara. Continuous-time non-symmetric algebraic Riccati theory: A matrix pencil approach. European Journal of Control, 18:74-81, 2012.

[4] D. Kremer and R. Stefan. Non-symmetric Riccati theory and linear quadratic Nash games. Notre Dame, MTNS 2002

nared(A1, A2, B1, B2, R, Q, S1 = 0, S2 = 0; scaling = 'B', pow2 = false, as = false, rtol::Real = nϵ, nrm = 1, indeig) -> (X, EVALS, F, Z, scalinfo, ind, HEVALS)

Computes X, the stabilizing solution (if as = false) or anti-stabilizing solution (if as = true) of the discrete-time nonsymmetric/nonhermitian algebraic Riccati equation

 A2'XA1 - X - (A2'XB1+S1)(R+B2'XB1)^(-1)(B2'XA1+S2') + Q = 0,

where A1 is an n×n matrix, A2 is an l×l matrix, and B1, B2, R, Q, S1, S2 are n×m, l×m, m×m, l×n, l×m, n×m matrices, respectively, with R nonsingular. Scalar-valued R is interpreted as an appropriately sized uniform scaling operator R*I. S1 and S2, if not specified, are set to zero matrices of appropriate sizes. The generalized Schur method of [1] and [2] is used. Solvents corrsponding to arbitrary collections of n generalized finite eigenvalues of the pair (H,J) (see below) can be selected by specifying their indices via indeig = ind, where ind must be an integer vector with at least n uniques indices being part of a permutation of the vector 1:n+l. For real data, the indices must additionally indicate a symmetric set of complex conjugate eigenvalues.

To enhance the accuracy of computations, a block oriented scaling of matrices R, Q and S1 is performed using the default setting scaling = 'B'. This scaling is performed only if norm(Q) > norm(B1)^2/norm(R). A general, eigenvalue computation oriented scaling combined with a block scaling is used if scaling = 'G' is selected. If A1 and A2 have the same sizes, an alternative, structure preserving scaling can be performed using the option scaling = 'S'. For the same case, a symmetric matrix equilibration based scaling can be employed if scaling = 'K', for which the underlying vector norm can be specified using the keyword argument nrm = p, where p = 1 is the default setting. Experimental structure preserving scalings can be performed using the options scaling = 'D', scaling = 'R' and scaling = 'T' (the last two only if A1 and A2 have the same sizes). Scaling can be disabled with the choice scaling = 'N'.

By default, the lower bound for the 1-norm reciprocal condition number rtol is n*ϵ, where n is the order of A1 and ϵ is the machine epsilon of the element type of A1.

X is a stabilizing or anti-stabilizing solvent of the Riccati equation for the selected/chosen finite generalized eigenvalues of the pair (H,J), where H = [ A1 0 B1; -Q I -S1; S2' 0 R ] and J = [ I 0 0; 0 A2' 0; 0 -B2' 0 ]. X = nothing if no such a solution exists, in which case Z, ind and HEVALS (see below) contain additional information to serve for possible alternative choices.

EVALS is a vector containing n stable or n anti-stable finite generalized eigenvalues of the pair (H,J), which underlies the computation of the solution. These are also the eigenvalues of A1-B1F.

F is the stabilizing or anti-stabilizing gain matrix F = (R+B2'XB1)^(-1)(B2'X*A1+S2').

Z = [U; V; W] is an orthogonal basis for the relevant stable/anti-stable deflating subspace such that X = Sx*(V/U)*Sxi and F = -Sr*(W/U)*Sxi, where Sx, Sxi and Sr are diagonal scaling matrices contained in the named tuple scalinfo as scalinfo.Sx, scalinfo.Sxi and scalinfo.Sr, respectively.

ind contains the indices of the eigenvalues used to compute the deflating subspace underlying the computation of the resulting solvent X.

HEVALS is a vector containing the n+l finite generalized eigenvalues of the pair (H,J). Only HEVALS can be computed with the call HEVALS = nared(A1, A2, B1, B2, R, Q, S1 = 0, S2 = 0; onlyeigvals = true, ...).

Note: The considered nonsymmetric/nonhermitian Riccati equation is related to the computation of Nash equilibria in non-cooperative differential games [2]. The option to select eigenvalues via indeig = ind is instrumental for determining various solvents.

Reference:

[1] W.F. Arnold, III and A.J. Laub, Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations, Proc. IEEE, 72:1746-1754, 1984.

[3] M. Jungers, C. Oara, H.-A. Khandil and R. Stefan. General matrix pencil techniques for solving discrete-time nonsymmetric algebraic Riccati equations, SIAM Journal on Matrix Analysis and Applications, 31:1257-1278, 2009.

garec(A, E, G, Q = 0; scaling = 'B', pow2 = false, as = false, rtol::Real = nϵ, nrm = 1) -> (X, EVALS, Z, scalinfo)

Compute X, the hermitian/symmetric stabilizing solution (if as = false) or anti-stabilizing solution (if as = true) of the generalized continuous-time algebraic Riccati equation

A'XE + E'XA - E'XGXE + Q = 0,

where G and Q are hermitian/symmetric matrices or uniform scaling operators and E is a nonsingular matrix. Scalar-valued G and Q are interpreted as appropriately sized uniform scaling operators G*I and Q*I. The generalized Schur method of [1] is used.

To enhance the accuracy of computations, a block scaling of matrices G and Q is performed, if the default setting scaling = 'B' is used. This scaling is however performed only if norm(Q) > norm(G). A general, eigenvalue computation oriented scaling combined with a block scaling is used if scaling = 'G' is selected. An alternative, structure preserving scaling can be performed using the option scaling = 'S'. A symmetric matrix equilibration based scaling is employed if scaling = 'K', for which the underlying vector norm can be specified using the keyword argument nrm = p, where p = 1 is the default setting. Scaling can be disabled with the choice scaling = 'N'. If pow2 = true, the scaling elements are enforced to the nearest power of 2 (default: pow2 = false).

By default, the lower bound for the 1-norm reciprocal condition number rtol is n*ϵ, where n is the order of A and ϵ is the machine epsilon of the element type of A.

EVALS is a vector containing the (stable or anti-stable) generalized eigenvalues of the pair (A-GXE,E).

Z = [U; V] is an orthogonal basis for the stable/anti-stable deflating subspace such that X = (Sx*(V/U)*Sxi)/E, where Sx and Sxi are diagonal scaling matrices contained in the named tuple scalinfo as scalinfo.Sx and scalinfo.Sxi, respectively.

Note: To solve the continuous-time algebraic Riccati equation

 A'XE + E'XA - E'XBR^(-1)B'XE + Q = 0,

with R a hermitian/symmetric matrix and B a compatible size matrix, G = BR^(-1)B' must be provided. This approach is not numerically suited when R is ill-conditioned and/or B has large norm.

Reference:

[1] W.F. Arnold, III and A.J. Laub, Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations, Proc. IEEE, 72:1746-1754, 1984.

garec(A, E, B, R, Q, S; scaling = 'B', pow2 = false, as = false, rtol::Real = nϵ, nrm = 1) -> (X, EVALS, F, Z, scalinfo)

Compute X, the hermitian/symmetric stabilizing solution (if as = false) or anti-stabilizing solution (if as = true) of the generalized continuous-time algebraic Riccati equation

A'XE + E'XA - (E'XB+S)R^(-1)(B'XE+S') + Q = 0,

where R and Q are hermitian/symmetric matrices such that R is nonsingular, and E is a nonsingular matrix. Scalar-valued R and Q are interpreted as appropriately sized uniform scaling operators R*I and Q*I. S, if not specified, is set to S = zeros(size(B)). The generalized Schur method of [1] is used.

To enhance the accuracy of computations, a block oriented scaling of matrices R, Q and S is performed using the default setting scaling = 'B'. This scaling is performed only if norm(Q) > norm(B)^2/norm(R). A general, eigenvalue computation oriented scaling combined with a block scaling is used if scaling = 'G' is selected. An alternative, structure preserving scaling can be performed using the option scaling = 'S'. A symmetric matrix equilibration based scaling is employed if scaling = 'K', for which the underlying vector norm can be specified using the keyword argument nrm = p, where p = 1 is the default setting. Experimental structure preserving scalings can be performed using the options scaling = 'D' or scaling = 'T'. Scaling can be disabled with the choice scaling = 'N'. If pow2 = true, the scaling elements are enforced to the nearest power of 2 (default: pow2 = false).

By default, the lower bound for the 1-norm reciprocal condition number rtol is n*ϵ, where n is the order of A and ϵ is the machine epsilon of the element type of A.

EVALS is a vector containing the (stable or anti-stable) generalized eigenvalues of the pair (A-BF,E).

F is the stabilizing/anti-stabilizing gain matrix F = R^(-1)(B'XE+S').

Z = [U; V; W] is an orthogonal basis for the relevant stable/anti-stable deflating subspace such that X = (Sx*(V/U)*Sxi)/E and F = -Sr*(W/U)*Sxi, where Sx, Sxi and Sr are diagonal scaling matrices contained in the named tuple scalinfo as scalinfo.Sx, scalinfo.Sxi and scalinfo.Sr, respectively.

Reference:

[1] W.F. Arnold, III and A.J. Laub, Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations, Proc. IEEE, 72:1746-1754, 1984.

Example

julia> using LinearAlgebra

julia> A = [-6. -2. 1.; 5. 1. -1; -4. -2. -1.]
3×3 Array{Float64,2}:
 -6.0  -2.0   1.0
  5.0   1.0  -1.0
 -4.0  -2.0  -1.0

julia> E = [10. 3. 0.; 0. 5. -1.; 0. 0. 10.]
3×3 Array{Float64,2}:
 10.0  3.0   0.0
  0.0  5.0  -1.0
  0.0  0.0  10.0

julia> B = [1. 2.; 2. 0.; 0. 1.]
3×2 Array{Float64,2}:
 1.0  2.0
 2.0  0.0
 0.0  1.0

julia> R = [1. 0.; 0. 5.]
2×2 Array{Float64,2}:
 1.0  0.0
 0.0  5.0

julia> X, CLSEIG, F = garec(A,E,B,R,2I);

julia> X
3×3 Array{Float64,2}:
  0.0502214   0.0284089   -0.0303703
  0.0284089   0.111219    -0.00259162
 -0.0303703  -0.00259162   0.0618395

julia> A'*X*E+E'*X*A-E'*X*B*inv(R)*B'*X*E+2I
3×3 Array{Float64,2}:
  1.55431e-15  -1.9984e-15   -3.33067e-15
 -1.77636e-15   1.33227e-15  -3.33067e-15
 -2.88658e-15  -3.21965e-15   1.11022e-15

julia> CLSEIG
3-element Array{Complex{Float64},1}:
  -0.6184265391601464 + 0.2913286844595737im
  -0.6184265391601464 - 0.2913286844595737im
 -0.21613059964451786 + 0.0im

julia> eigvals(A-B*F,E)
3-element Array{Complex{Float64},1}:
 -0.6184265391601462 - 0.29132868445957383im
 -0.6184265391601462 + 0.2913286844595739im
  -0.216130599644518 + 0.0im

garec(A, E, B, G, R, Q, S; scaling = 'B', pw2 = false, as = false, rtol::Real = nϵ, nrm = 1) -> (X, EVALS, F, Z, scalinfo)

Compute X, the hermitian/symmetric stabilizing solution (if as = false) or anti-stabilizing solution (if as = true) of the generalized continuous-time algebraic Riccati equation

A'XE + E'XA - E'XGXE - (E'XB+S)R^(-1)(B'XE+S') + Q = 0,

where G, Q and R are hermitian/symmetric matrices such that R is nonsingular, and E is a nonsingular matrix. Scalar-valued G, R and Q are interpreted as appropriately sized uniform scaling operators G*I, R*I and Q*I. The generalized Schur method of [1] is used.

To enhance the accuracy of computations, a block oriented scaling of matrices G, R, Q and S is performed using the default setting scaling = 'B'. This scaling is performed only if norm(Q) > max(norm(G), norm(B)^2/norm(R)). A general, eigenvalue computation oriented scaling combined with a block scaling is used if scaling = 'G' is selected. An alternative, structure preserving scaling can be performed using the option scaling = 'S'. A symmetric matrix equilibration based scaling is employed if scaling = 'K', for which the underlying vector norm can be specified using the keyword argument nrm = p, where p = 1 is the default setting. Experimental structure preserving scalings can be performed using the options scaling = 'D' or scaling = 'T'. Scaling can be disabled with the choice scaling = 'N'. If pow2 = true, the scaling elements are enforced to the nearest power of 2 (default: pow2 = false).

By default, the lower bound for the 1-norm reciprocal condition number rtol is n*ϵ, where n is the order of A and ϵ is the machine epsilon of the element type of A.

EVALS is a vector containing the (stable or anti-stable) generalized eigenvalues of the pair (A-BF-GXE,E).

F is the stabilizing/anti-stabilizing gain matrix F = R^(-1)(B'XE+S').

Z = [U; V; W] is an orthogonal basis for the relevant stable/anti-stable deflating subspace such that X = (Sx*(V/U)*Sxi)/E and F = -Sr*(W/U)*Sxi, where Sx, Sxi and Sr are diagonal scaling matrices contained in the named tuple scalinfo as scalinfo.Sx, scalinfo.Sxi and scalinfo.Sr, respectively.

Reference:

[1] W.F. Arnold, III and A.J. Laub, Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations, Proc. IEEE, 72:1746-1754, 1984.

gared(A, E, B, R, Q, S; scaling = 'B', pow2 = false, as = false, rtol::Real = nϵ, nrm = 1) -> (X, EVALS, F, Z, scalinfo)

Compute X, the hermitian/symmetric stabilizing solution (if as = false) or anti-stabilizing solution (if as = true) of the generalized discrete-time algebraic Riccati equation

A'XA - E'XE - (A'XB+S)(R+B'XB)^(-1)(B'XA+S') + Q = 0,

where R and Q are hermitian/symmetric matrices, and E ist non-singular. Scalar-valued R and Q are interpreted as appropriately sized uniform scaling operators R*I and Q*I. S, if not specified, is set to S = zeros(size(B)).

To enhance the accuracy of computations, a block oriented scaling of matrices R, Q and S is performed using the default setting scaling = 'B'. This scaling is performed only if norm(Q) > norm(B)^2/norm(R). A general, eigenvalue computation oriented scaling combined with a block scaling is used if scaling = 'G' is selected. An alternative, structure preserving scaling can be performed using the option scaling = 'S'. A symmetric matrix equilibration based scaling is employed if scaling = 'K', for which the underlying vector norm can be specified using the keyword argument nrm = p, where p = 1 is the default setting. Experimental structure preserving scalings can be performed using the options scaling = 'D', scaling = 'R' and scaling = 'T'. Scaling can be disabled with the choice scaling = 'N'. If pow2 = true, the scaling elements are enforced to the nearest power of 2 (default: pow2 = false).

By default, the lower bound for the 1-norm reciprocal condition number rtol is n*ϵ, where n is the order of A and ϵ is the machine epsilon of the element type of A.

EVALS is a vector containing the (stable or anti-stable) generalized eigenvalues of the pair (A-BF,E).

F is the stabilizing/anti-stabilizing gain matrix F = (R+B'XB)^(-1)(B'XA+S').

Z = [U; V; W] is an orthogonal basis for the relevant stable/anti-stable deflating subspace such that X = (Sx*(V/U)*Sxi)/E and F = -Sr*(W/U)*Sxi, where Sx, Sxi and Sr are diagonal scaling matrices contained in the named tuple scalinfo as scalinfo.Sx, scalinfo.Sxi and scalinfo.Sr, respectively.

Reference:

[1] W.F. Arnold, III and A.J. Laub, Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations, Proc. IEEE, 72:1746-1754, 1984.

Example

julia> using LinearAlgebra

julia> A = [-6. -2. 1.; 5. 1. -1; -4. -2. -1.]
3×3 Array{Float64,2}:
 -6.0  -2.0   1.0
  5.0   1.0  -1.0
 -4.0  -2.0  -1.0

julia> E = [10. 3. 0.; 0. 5. -1.; 0. 0. 10.]
3×3 Array{Float64,2}:
 10.0  3.0   0.0
  0.0  5.0  -1.0
  0.0  0.0  10.0

julia> B = [1. 2.; 2. 0.; 0. 1.]
3×2 Array{Float64,2}:
 1.0  2.0
 2.0  0.0
 0.0  1.0

julia> R = [1. 0.; 0. 5.]
2×2 Array{Float64,2}:
 1.0  0.0
 0.0  5.0

julia> X, CLSEIG, F = gared(A,E,B,R,2I);

julia> X
3×3 Array{Float64,2}:
  0.065865   -0.0147205  -0.0100407
 -0.0147205   0.0885939   0.0101422
 -0.0100407   0.0101422   0.0234425

julia> A'*X*A-E'*X*E-A'*X*B*inv(R+B'*X*B)*B'*X*A+2I
3×3 Array{Float64,2}:
 -1.33227e-15  -2.48412e-15   1.38778e-16
 -2.498e-15    -4.44089e-16  -6.50521e-16
  1.80411e-16  -5.91541e-16  -1.33227e-15

julia> CLSEIG
3-element Array{Complex{Float64},1}:
  -0.084235615751339 - 0.0im
  -0.190533552034239 - 0.0im
 -0.5238922629921539 - 0.0im

julia> eigvals(A-B*F,E)
3-element Array{Float64,1}:
 -0.5238922629921539
 -0.19053355203423886
 -0.08423561575133902

uqme(A, B, C; as = true, disc = true, scaling = 'B', pow2 = false,  
     rtol::Real = nϵ, nrm = 1, indeig) -> (X, EVALS, Z, scalinfo, ind, HEVALS)

Compute a solvent X of the unilateral quadratic matrix equation (UQMU)

 AX^2 + BX + C = 0,

where A, B and C are n×n matrices, and X is an n×n matrix. The generalized Schur-form based technique of [1] is used, where deflating subspaces [U; V] corresponding to the selected n finite generalized eigenvalues of a suitable pencil F-λG are determined and X = V*inv(U) is the resulting solvent. A preliminary deflating of infinite eigenvalues is performed in the case of a singular A. If as = true and disc = true the dominant solvent is attempted to be computed (default), corresponding to the eigenvalues with the largest moduli, while if as = false and disc = true, the minimal solvent is targetted, corresponding to the eigenvalues with the smallest moduli. Alternative solvents can be targetted by selecting eigenvalues with largest real parts with as = true and disc = false or smallest real parts using as = false and disc = false. Solvents corresponding to arbitrary collections of eigenvalues can be selected by specifying their indices via indeig = ind, where ind must be an integer vector with at least n uniques indices being part of a permutation of the vector 1:nf, where nf is the number of finite eigenvalues. For real data, the indices must additionally indicate a symmetric set of complex conjugate eigenvalues. The case A = I is explicitly handled, by solving a (simpler) nonsymmetric/nonhermitian algebraic Riccati equation.

To enhance the accuracy of computations, a block scaling of matrices I and C is performed, if the default setting scaling = 'B' is used. This scaling is however performed only if norm(C) > 1. A general, eigenvalue computation oriented scaling combined with a block scaling is used if scaling = 'G' is selected. An alternative, experimental structure preserving scaling can be performed using the option scaling = 'S'. A symmetric matrix equilibration based scaling is employed if scaling = 'K', for which the underlying vector norm can be specified using the keyword argument nrm = p, where p = 1 is the default setting. Scaling can be disabled with the choice scaling = 'N'. If pow2 = true, the scaling elements are enforced to the nearest power of 2 (default: pow2 = false).

By default, the lower bound for the 1-norm reciprocal condition number rtol is n*ϵ, where n is the order of A and ϵ is the machine epsilon of the element type of C. For rank determinations, rank revealing QR-decomposition with column pivoting is used, where besides rtol, and absolute tolerance atol (Default: atol = 0), is also used.

X is a solvent of the quadratic matrix equation for the selected/chosen finite generalized eigenvalues of the pair (F,G), where F = [0 I; -C -B] and G = [I 0; 0 A]. X = nothing if no such a solution exists, in which case Z, ind and HEVALS (see below) contain additional information to serve for possible alternative choices.

EVALS is a vector containing n generalized eigenvalues of the pair (F,G), which underlies the computation of the solution. These are also the eigenvalues of X. The selection of the eigenvalues is performed as follows. If disc = true, then the n eigenvalues with the smallest moduli are selected if as = false, or with the largest moduli if as = true. If disc = false, then the n eigenvalues with the smallest real parts are selected if as = false, or with the largest real parts if as = true. Possible permutations of eigenvalues can occur if the n-th eigenvalue belongs to a complex conjugated pair. The infinite eigenvalues (e.g., if A is singular) are excluded from any selection.

Z = [U; V] is an orthogonal basis for thenarec(A1, A2, B1, B2, R, Q, S1 = 0, S2 = 0; deflating subspace corresponding to EVALS, such that X = Sx*(V/U)*Sxi, where Sx and Sxi are diagonal scaling matrices contained in the named tuple scalinfo as scalinfo.Sx and scalinfo.Sxi, respectively.

ind contains the indices of the finite eigenvalues used to compute the deflating subspace underlying the computation of the resulting solvent X.

HEVALS is a vector containing the 2n generalized eigenvalues of the pair (F,G). Only HEVALS can be computed with the call HEVALS = uqme(A, B, C; onlyeigvals = true, ...).

Note: Generically, for complex data, there are up most C(2n,n) complex solvents which correspond to all possible combinations of n eigenvalues of the 2n eigenvalues of the pencil F-λG. For real data, there are usually less solvents than for complex data. All solvents can be determined by generating an exhaustive set of C(2n,n) permutations and checking for each case if a solution exists (see some of the attached tests examples using try/catch constructs).

Reference:

[1] N. J. Higham and H.-M. Kim, Numerical analysis of a quadratic matrix equation. IMA Journal of Numerical Analysis 20:499–519, 2000.