andreasvarga / MatrixEquations.jl

Solution of Lyapunov, Sylvester and Riccati matrix equations using Julia
MIT License
82 stars 8 forks source link

Condition Matrix(T)' = Matrix(T') not fulfilled for all defined operators #2

Closed andreasvarga closed 3 years ago

andreasvarga commented 4 years ago

Let T be, for example, a Lyapunov operator T:X -> Y, which maps the upper triangular parts of symmetric matrices X to the upper triangular parts of the symmetric matrices Y computed as Y := AX+XA'. Similarly, the transposed operator Tt =T' maps X to Y := A'X+XA. Consider the example:

julia> A = [11 12; 21 22]
2×2 Array{Int64,2}:
 11  12
 21  22

julia> T = lyapop(A,her=true)
Linear operator
  nrow: 3
  ncol: 3
  eltype: Int64
  symmetric: false
  hermitian: false

julia> Tt = lyapop(A',her=true)
Linear operator
  nrow: 3
  ncol: 3
  eltype: Int64
  symmetric: false
  hermitian: false

julia> Matrix(T)
3×3 Array{Int64,2}:
 22  24   0
 21  33  12
  0  42  44

julia> Matrix(T')
3×3 Array{Int64,2}:
 22  42   0
 12  33  21
  0  24  44

julia> Matrix(Tt)
3×3 Array{Int64,2}:
 22  42   0
 12  33  21
  0  24  44

julia> Matrix(T)' == Matrix(T')
false

julia> Matrix(Tt)' == Matrix(Tt')
false

However we also have, as expected,

julia> Matrix(T') == Matrix(Tt)
true

julia> Matrix(Tt') == Matrix(T)
true

Fixing this inconsistent behaviour would be very desirable in a future release.

dpo commented 4 years ago

Are you defining Matrix(op) where op is a Lyapunov operator anywhere? Or are you relying on the default implementation of Matrix()?

andreasvarga commented 4 years ago

I am using the default implementation.

Dominique notifications@github.com schrieb am Di., 3. Dez. 2019, 01:57:

Are you defining Matrix(op) where op is a Lyapunov operator anywhere? Or are you relying on the default implementation of Matrix()?

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/andreasvarga/MatrixEquations.jl/issues/2?email_source=notifications&email_token=ALJDHEGADYAIMPIYBNZGRPDQWWVHRA5CNFSM4JOFBIY2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEFXWOAY#issuecomment-560948995, or unsubscribe https://github.com/notifications/unsubscribe-auth/ALJDHEB3CO45BSONTZ3F4PTQWWVHRANCNFSM4JOFBIYQ .

andreasvarga commented 3 years ago

The Matrix(T)' == Matrix(T') issue has been solved by determining Matrix(T') as Matrix(T)' (not a very elegant solution I admit!)