Riccati Matrix Equation Solvers
Standard Riccati Matrix Equations
MatrixEquations.arec — Functionarec(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.0imarec(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.0imarec(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.
MatrixEquations.ared — Functionared(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.5Generalized Riccati Matrix Equations
MatrixEquations.garec — Functiongarec(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.0imgarec(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.
MatrixEquations.gared — Functiongared(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