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)

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. 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 A1.

EVALS is a vector containing n stable or n anti-stable eigenvalues of the Hamiltonian-like matrix H := [A1 -G; -Q -A2], 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.

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.

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)

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 generalized eigenvalues of a suitable pencil F-λG are determined and X = V*inv(U) is the resulting solvent. 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 corrsponding 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:2*n. 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.

EVALS is a vector containing n generalized eigenvalues of the pencil F-λG, where F = [0 I; -C -B] and G = [I 0; 0 A]. 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 exact infinite eigenvalues (e.g., if A is singular) are excluded from any selection.

Z = [U; V] is an orthogonal basis for the 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 eigenvalues used to compute the deflating subspace underlying the computation of the resulting solvent X.

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.