Sylvester Matrix Equation Solvers

Standard Sylvester Matrix Equations

MatrixEquations.sylvcFunction
X = sylvc(A,B,C)

Solve the continuous Sylvester matrix equation

            AX + XB = C

using the Bartels-Stewart Schur form based approach. A and B are square matrices, and A and -B must not have common eigenvalues.

The following particular cases are also adressed:

X = sylvc(α*I,B,C)  or  X = sylvc(α,B,C)

Solve the matrix equation X(αI+B) = C.

X = sylvc(A,β*I,C)  or  X = sylvc(A,β,C)

Solve the matrix equation (A+βI)X = C.

X = sylvc(α*I,β*I,C)  or  sylvc(α,β,C)

Solve the matrix equation (α+β)X = C.

x = sylvc(α,β,γ)

Solve the equation (α+β)x = γ.

Example

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

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

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

julia> X = sylvc(A, B, C)
2×2 Array{Float64,2}:
 -4.46667   1.93333
  3.73333  -1.8

julia> A*X + X*B - C
2×2 Array{Float64,2}:
  2.66454e-15  1.77636e-15
 -3.77476e-15  4.44089e-16
MatrixEquations.sylvcs!Function
sylvcs!(A,B,C; adjA = false, adjB = false)

Solve the continuous Sylvester matrix equation

            op(A)X + Xop(B) =  C,

where op(A) = A or op(A) = A' if adjA = false or adjA = true, respectively, and op(B) = B or op(B) = B' if adjB = false or adjB = true, respectively. A and B are square matrices in Schur forms, and A and -B must not have common eigenvalues. C contains on output the solution X.

MatrixEquations.sylvdFunction
X = sylvd(A,B,C)

Solve the discrete Sylvester matrix equation

            AXB + X = C

using an extension of the Bartels-Stewart Schur form based approach. A and B are square matrices, and A and -B must not have common reciprocal eigenvalues.

The following particular cases are also adressed:

X = sylvd(α*I,B,C)  or  X = sylvd(α,B,C)

Solve the matrix equation X(αB+I) = C.

X = sylvd(A,β*I,C)   or  X = sylvd(A,β,C)

Solve the matrix equation (βA+I)X = C.

X = sylvd(α*I,β*I,C)  or  X = sylvd(α,β,C)

Solve the matrix equation (αβ+1)X = C.

x = sylvd(α,β,γ)

Solve the equation (αβ+1)x = γ.

Example

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

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

julia> C = [-1. -2.; 2. -1.]
2×2 Array{Float64,2}:
 -1.0  -2.0
  2.0  -1.0
   
julia> X = sylvd(A, B, C)
2×2 Array{Float64,2}:
 -2.46667  -2.73333
  2.4       1.86667

julia> A*X*B + X - C
2×2 Array{Float64,2}:
  8.88178e-16   8.88178e-16
 -3.9968e-15   -5.55112e-15
MatrixEquations.sylvds!Function
sylvds!(A,B,C; adjA = false, adjB = false)

Solve the discrete Sylvester matrix equation

            op(A)Xop(B) + X =  C,

where op(A) = A or op(A) = A' if adjA = false or adjA = true, respectively, and op(B) = B or op(B) = B' if adjB = false or adjB = true, respectively. A and B are square matrices in Schur forms, and A and -B must not have common reciprocal eigenvalues. C contains on output the solution X.

Generalized Sylvester Matrix Equations

MatrixEquations.gsylvFunction
X = gsylv(A,B,C,D,E)

Solve the generalized Sylvester matrix equation

          AXB + CXD = E

using a generalized Schur form based approach. A, B, C and D are square matrices. The pencils A-λC and D+λB must be regular and must not have common eigenvalues.

The following particular cases are also adressed:

X = gsylv(A,B,E)

Solve the generalized Sylvester matrix equation AXB = E.

X = gsylv(A,B,γ*I,E)  or  X = gsylv(A,B,γ,E)

Solve the generalized Sylvester matrix equation AXB +γX = E.

X = gsylv(A,B,γ*I,D,E)  or  X = gsylv(A,B,γ,D,E)

Solve the generalized Sylvester matrix equation AXB +γXD = E.

X = gsylv(A,B,C,δ*I,E)  or  X = gsylv(A,B,C,δ,E)

Solve the generalized Sylvester matrix equation AXB +CXδ = E.

Example

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

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

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

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

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

julia> X = gsylv(A, B, C, D, E)
2×2 Array{Float64,2}:
 -0.52094   -0.0275792
 -0.168539   0.314607
 
julia> A*X*B + C*X*D - E
2×2 Array{Float64,2}:
 4.44089e-16  8.88178e-16
 6.66134e-16  0.0
MatrixEquations.gsylvs!Function
X = gsylvs!(A,B,C,D,E; adjAC=false, adjBD=false, CASchur = false, DBSchur = false)

Solve the generalized Sylvester matrix equation

            op1(A)Xop2(B) + op1(C)Xop2(D) = E,

where A, B, C and D are square matrices, and

op1(A) = A and op1(C) = C if adjAC = false;

op1(A) = A' and op1(C) = C' if adjAC = true;

op2(B) = B and op2(D) = D if adjBD = false;

op2(B) = B' and op2(D) = D' if adjBD = true.

The matrix pair (A,C) is in a generalized real or complex Schur form. The matrix pair (B,D) is in a generalized real or complex Schur form if DBSchur = false or the matrix pair (D,B) is in a generalized real or complex Schur form if DBSchur = true. The pencils A-λC and D+λB must be regular and must not have common eigenvalues.

Sylvester Systems of Matrix Equation

MatrixEquations.sylvsysFunction
(X,Y) = sylvsys(A,B,C,D,E,F)

Solve the Sylvester system of matrix equations

            AX + YB = C
            DX + YE = F,

where (A,D), (B,E) are pairs of square matrices of the same size. The pencils A-λD and -B+λE must be regular and must not have common eigenvalues.

Example

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

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

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

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

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

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

julia> X, Y = sylvsys(A, B, C, D, E, F);

julia> X
2×2 Array{Float64,2}:
  1.388  -1.388
 -0.892   0.892

julia> Y
2×2 Array{Float64,2}:
 -1.788  0.192
  0.236  0.176

julia> A*X + Y*B - C
2×2 Array{Float64,2}:
  6.66134e-16  2.22045e-15
 -3.10862e-15  2.66454e-15

julia> D*X + Y*E - F
2×2 Array{Float64,2}:
  1.33227e-15  -2.22045e-15
 -4.44089e-16   4.44089e-16
MatrixEquations.sylvsyss!Function
(X,Y) = sylvsyss!(A,B,C,D,E,F)

Solve the Sylvester system of matrix equations

            AX + YB = C
            DX + YE = F,

where (A,D), (B,E) are pairs of square matrices of the same size in generalized Schur forms. The pencils A-λD and -B+λE must be regular and must not have common eigenvalues. The computed solution (X,Y) is contained in (C,F).

Note: This is an enhanced interface to the LAPACK.tgsyl! function to also cover the case when A, B, D and E are real matrices and C and F are complex matrices.

MatrixEquations.dsylvsysFunction
(X,Y) = dsylvsys(A,B,C,D,E,F)

Solve the dual Sylvester system of matrix equations

   AX + DY = C
   XB + YE = F ,

where (A,D), (B,E) are pairs of square matrices of the same size. The pencils A-λD and -B+λE must be regular and must not have common eigenvalues.

Example

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

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

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

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

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

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

julia> X, Y = dsylvsys(A, B, C, D, E, F);

julia> X
2×2 Array{Float64,2}:
  2.472  -1.648
 -1.848   1.232

julia> Y
2×2 Array{Float64,2}:
 -0.496  -0.336
  0.264   0.824

julia> A*X + D*Y - C
2×2 Array{Float64,2}:
  4.44089e-16  0.0
 -3.55271e-15  1.55431e-15
 
julia> X*B + Y*E - F
2×2 Array{Float64,2}:
 -8.88178e-16   0.0
  8.88178e-16  -4.44089e-16
MatrixEquations.dsylvsyss!Function
(X,Y) = dsylvsyss!(A,B,C,D,E,F)

Solve the dual Sylvester system of matrix equations

A'X + D'Y = C
XB' + YE' = F,

where (A,D), (B,E) are pairs of square matrices of the same size in generalized Schur forms. The pencils A-λD and -B+λE must be regular and must not have common eigenvalues. The computed solution (X,Y) is contained in (C,F).