Problem with lyap(Q,R)

ylish · July 16, 2021, 12:36pm

Hello everybody,

I am new to Julia and I am trying to rewrite and improve some code that I have in MATLAB. I have come across something that I cannot fix. Could someone see what the issue is?

I have the following code:

using LinearAlgebra

A = [1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0]

B = [0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; -1252.0 1250.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0; 1250.0 -1252.0 2.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0; 0.0 2.0 -1252.0 1250.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0; 0.0 0.0 1250.0 -1252.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0; -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1252.0 -1250.0 0.0 0.0; 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1250.0 1252.0 -2.0 0.0; 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -2.0 1252.0 -1250.0; 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1250.0 1252.0; 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0]

lyap(Matrix(B’),A)

but I get the following message:
LAPACKException(1)

Stacktrace:
[1] chklapackerror(ret::Int64)
@ LinearAlgebra.LAPACK /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lapack.jl:38
[2] trsyl!(transa::Char, transb::Char, A::Matrix{Float64}, B::Matrix{Float64}, C::Matrix{Float64}, isgn::Int64)
@ LinearAlgebra.LAPACK /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lapack.jl:6463
[3] trsyl!
@ /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lapack.jl:6444 [inlined]
[4] lyap(A::Matrix{Float64}, C::Matrix{Float64})
@ LinearAlgebra /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/dense.jl:1567
[5] top-level scope
@ In[17]:4
[6] eval
@ ./boot.jl:360 [inlined]
[7] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:1094

In MATLAB there is no issue with these two matrices and the same exact comand. Would someone know what the issue could be ? ][=

Sukera · July 16, 2021, 12:55pm

Welcome! In order to make it easier to help you, please wrap your code in triple backtickes (```) to format it as code. For more information, please read Make it easier to help you. Please also make sure you’re aware of the Differences from MATLAB.

I’m not sure why this works in MATLAB in the first place - B' has negative complex conjugate eigen values, which are forbidden per the dosctring of lyap (or I’m misunderstanding something about the dosctring):

julia> eigvals(B')                         
16-element Vector{ComplexF64}:             
 -0.6761400549933252 - 0.978143841175781im 
 -0.6761400549933252 + 0.978143841175781im 
 -0.5754555540509579 - 1.6824830146993683im
 -0.5754555540509579 + 1.6824830146993683im
 -0.5000999500060193 - 50.007504436433635im
 -0.5000999500060193 + 50.007504436433635im
  -0.500099870149834 - 50.027497437709776im
  -0.500099870149834 + 50.027497437709776im
  0.5000998701498371 - 50.027497437709776im
  0.5000998701498371 + 50.027497437709776im
  0.5000999500060104 - 50.00750443643361im 
  0.5000999500060104 + 50.00750443643361im 
  0.5754555540509454 - 1.6824830146993628im
  0.5754555540509454 + 1.6824830146993628im
   0.676140054993325 - 0.9781438411758407im
   0.676140054993325 + 0.9781438411758407im

help?> lyap                                                                                                                                                                                 
search: lyap displayable                                                                                                                                                                    
                                                                                                                                                                                            
  lyap(A, C)                                                                                                                                                                                
                                                                                                                                                                                            
  Computes the solution X to the continuous Lyapunov equation AX + XA' + C = 0,
  where no eigenvalue of A has a zero real part and no two eigenvalues are
  negative complex conjugates of each other.

zdenek_hurak · July 16, 2021, 1:24pm

You may also want to try some alternative solvers. First, let me - just for convenience - rename your matrices so that they correspond to the common conventions used by solvers. Your A will become C, your B will become A after transposition:

C = [1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0]
A = [0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; -1252.0 1250.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0; 1250.0 -1252.0 2.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0; 0.0 2.0 -1252.0 1250.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0; 0.0 0.0 1250.0 -1252.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0; -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1252.0 -1250.0 0.0 0.0; 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1250.0 1252.0 -2.0 0.0; 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -2.0 1252.0 -1250.0; 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1250.0 1252.0; 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0]'

Now, I am aware of two alternatives for solving Lyapunov equation:

julia> using MatrixEquations

julia> X = lyapc(A,C)
16×16 Matrix{Float64}:
  15279.3          -15293.4           12471.6          …  -223.391        -119.321         115.15
 -15293.4           15299.2          -12479.8              227.77          114.793        -119.321
  12471.6          -12479.8           15299.2              114.793         227.77         -223.391
 -12455.4           12471.6          -15293.4             -119.321        -223.391         227.961
     -2.78533e-12       2.2311e-11       -2.34454e-8         1.21925e-11    -9.91488e-10     9.90487e-10
      1.06155e-11      -2.97575e-11       2.34695e-8   …     0.5             9.91008e-10    -9.90498e-10
      ⋮                                                ⋱                                     ⋮
      9.78416e-10      -9.78896e-10      -3.18567e-11        9.85301e-12     2.73115e-14     1.40998e-14
    227.961          -223.391          -119.321              7.10242        -6.05262         3.9942
   -223.391           227.77            114.793             -4.9365          4.00557        -6.05262
   -119.321           114.793           227.77               4.00557        -4.9365          7.10242
    115.15           -119.321          -223.391        …    -6.05262         7.10242        -4.92681

julia> A*X+X*A'+C
16×16 Matrix{Float64}:
  3.3527e-8    -3.49999e-8  -1.3459e-9   -9.35604e-10  …   2.41243e-9   -2.72769e-9    4.49009e-9   -2.63853e-9
 -3.49999e-8    3.64698e-8   1.17352e-9   1.10988e-9      -2.49065e-9    2.80954e-9   -4.42813e-9    2.57208e-9
 -1.3459e-9     1.17352e-9  -4.86263e-8   4.56729e-8       3.24606e-9   -3.6765e-9     4.3874e-9    -2.1422e-9
 -9.35604e-10   1.10988e-9   4.56729e-8  -4.27212e-8      -3.22284e-9    3.65187e-9   -4.42378e-9    2.18074e-9
 -2.49747e-9    2.48474e-9  -1.66074e-9   1.65346e-9      -7.63123e-11   8.03766e-11   3.14344e-11  -3.66214e-11
  5.33328e-9   -5.326e-9     4.77667e-9  -4.77121e-9   …   7.03722e-11  -7.77618e-11   3.45324e-12   5.57066e-12
  ⋮                                                    ⋱                                             ⋮
  1.88513e-9   -1.88631e-9   1.96059e-9  -1.95914e-9       2.40963e-11  -2.31539e-11  -1.49818e-11   1.37126e-11
  2.41243e-9   -2.49065e-9   3.24606e-9  -3.22284e-9      -7.75162e-11   8.55078e-11   1.25947e-10  -1.4731e-10
 -2.72769e-9    2.80954e-9  -3.6765e-9    3.65187e-9       8.55078e-11  -9.43823e-11  -1.29694e-10   1.51797e-10
  4.49009e-9   -4.42813e-9   4.3874e-9   -4.42378e-9       1.25947e-10  -1.29694e-10  -1.10189e-10   1.43357e-10
 -2.63853e-9    2.57208e-9  -2.1422e-9    2.18074e-9   …  -1.4731e-10    1.51797e-10   1.43357e-10  -1.77527e-10

and

julia> using ControlMatrixEquations

julia> X = lyapc(A,C)
16×16 Matrix{Float64}:  
15279.3          -15293.3           12471.6          …  -223.387        -119.323         115.147 -15293.3           15299.2          -12479.8              227.773         114.79         -119.323  12471.6          -12479.8           15299.2              114.79          227.773        -223.387 -12455.4           12471.6          -15293.3             -119.323        -223.387         227.964     -2.61548e-12       2.21302e-11      -2.34452e-8         1.21753e-11    -9.91495e-10     9.90495e-10      1.04541e-11      -2.9604e-11        2.34692e-8   …     0.5             9.91015e-10    -9.90506e-10      ⋮                                                ⋱                                     ⋮      9.78423e-10      -9.78905e-10      -3.18727e-11        9.853e-12       2.7605e-14      1.39602e-14
    227.964          -223.387          -119.323              7.10262        -6.0528          3.99402
   -223.387           227.773           114.79              -4.93629         4.00538        -6.0528
   -119.323           114.79            227.773              4.00538        -4.93629         7.10262
    115.147          -119.323          -223.387        …    -6.0528          7.10262        -4.92661

julia> A*X+X*A'+C16×16 Matrix{Float64}:  3.36389e-8   -3.51083e-8  -1.50165e-9  -8.36706e-10  …   2.40959e-9   -2.72439e-9    4.48801e-9   -2.64033e-9 -3.50854e-8    3.65519e-8   1.28263e-9   1.05764e-9      -2.48701e-9    2.80511e-9   -4.42472e-9    2.57256e-9 -1.57833e-9    1.30631e-9  -4.85399e-8   4.55455e-8       3.24075e-9   -3.67374e-9    4.38195e-9   -2.14322e-9 -7.61101e-10   1.03502e-9   4.56016e-8  -4.2609e-8       -3.21617e-9    3.64777e-9   -4.42071e-9    2.18392e-9 -2.49929e-9    2.49383e-9  -1.66438e-9   1.6571e-9       -7.63976e-11   8.0405e-11    3.1477e-11   -3.67208e-11  5.3351e-9    -5.32782e-9   4.7803e-9   -4.77485e-9   …   7.05427e-11  -7.78471e-11   3.31113e-12   5.65592e-12  ⋮                                                    ⋱                                             ⋮  1.88534e-9   -1.88663e-9   1.96076e-9  -1.95922e-9       2.41025e-11  -2.31735e-11  -1.49392e-11   1.36744e-11
  2.41596e-9   -2.49111e-9   3.24749e-9  -3.21498e-9      -7.74358e-11   8.5367e-11    1.25663e-10  -1.47177e-10
 -2.72917e-9    2.80772e-9  -3.67852e-9   3.64461e-9       8.53111e-11  -9.41252e-11  -1.29703e-10   1.51957e-10
  4.49071e-9   -4.42183e-9   4.38627e-9  -4.41719e-9       1.25676e-10  -1.2959e-10   -1.10448e-10   1.43313e-10
 -2.63916e-9    2.5658e-9   -2.14402e-9   2.17699e-9   …  -1.47204e-10   1.5186e-10    1.4371e-10   -1.77583e-10

fph · July 16, 2021, 1:30pm

Maybe I am missing something, but it doesn’t seem that it is the case, to me. “Negative complex conjugate eigenvalues”, in this context, should mean “there are no two eigenvalues such that \lambda_i + \overline{\lambda_j} = 0”. Some of those sums are very close to zero, but not exactly equal, it seems to me.

EDIT: I just checked the docs for dtrsyl.f, that’s what the error code means:

*          = 1: A and B have common or very close eigenvalues; perturbed
*               values were used to solve the equation (but the matrices
*               A and B are unchanged).

So (1) that INFO value appears even if the eigenvalues are merely close, not exactly equal, and (2) That INFO value alone should not cause a failure, since a valid solution is returned anyway (although this means that the problem is ill-conditioned and the returned solution is likely garbage).

Note that Matlab does not use trsyl from Lapack, but routines from Slicot, so it does not run the same code.

Sukera · July 16, 2021, 1:44pm

I see - in that case, seems like the function handling that “error” is misbehaving here:

From LinearAlgebra/src/lapack.jl:37

function chklapackerror(ret::BlasInt)                                   
    if ret == 0                                                         
        return                                                          
    elseif ret < 0                                                      
        throw(ArgumentError("invalid argument #$(-ret) to LAPACK call"))
    else # ret > 0                                                      
        throw(LAPACKException(ret))                                     
    end                                                                 
end

And it shouldn’t necessarily always throw here (or maybe a different error function should be used in LAPACK.trsyl!). Please open an issue on github about this!

rdeits · July 16, 2021, 2:21pm

What Julia version are you using? In Julia 1.6.1 your exact example appears to work correctly:

julia> lyap(Matrix(B'),A)
16×16 Matrix{Float64}:
  284.002        -283.488         234.109
...

julia> X = lyap(Matrix(B'), A);

julia> isapprox(B' * X + X * B + A, zeros(size(A)), atol=1e-8)
true

gbaraldi · July 16, 2021, 2:27pm

Thats weird, for me it fails on rosetta

julia> versioninfo()
Julia Version 1.6.1
Commit 6aaedecc44 (2021-04-23 05:59 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin18.7.0)
  CPU: Apple M1
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-11.0.1 (ORCJIT, westmere)
Environment:
  JULIA_NUM_THREADS = 4

And on native

versioninfo()
Julia Version 1.8.0-DEV.138
Commit 018977209b* (2021-07-06 11:28 UTC)
Platform Info:
  OS: macOS (arm64-apple-darwin20.5.0)
  CPU: Apple M1
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-12.0.0 (ORCJIT, cyclone)
Environment:
  JULIA_NUM_THREADS = 4

Sukera · July 16, 2021, 2:35pm

It definitely also fails for me on

julia> versioninfo()                           
Julia Version 1.7.0-DEV.1072                   
Commit 364660c5fd (2021-05-07 06:23 UTC)       
Platform Info:                                 
  OS: Linux (x86_64-linux-gnu)                 
  CPU: Intel(R) Core(TM) i7-6600U CPU @ 2.60GHz
  WORD_SIZE: 64                                
  LIBM: libopenlibm                            
  LLVM: libLLVM-11.0.1 (ORCJIT, skylake)       
Environment:                                   
  JULIA_PKG_SERVER =                           
  JULIA_NUM_THREADS = 4

ylish · July 16, 2021, 8:19pm

Thank you very much! I have marked this as the solution as it might help others as well looking to solve these equations.

ylish · July 16, 2021, 8:20pm

Thank you for the help! Were you asking me to open an issue on github? I haven’t done so before thus could you let me know what I should do?

ylish · July 16, 2021, 8:21pm

Thank you! This makes sense!

gbaraldi · July 16, 2021, 8:34pm

Not sure if you want to open this one but go here Issues · JuliaLang/julia · GitHub click new issue. Also give a read to julia/CONTRIBUTING.md at e1aeb8a95ede8212ea09bdb49bdd805f66dfb9d5 · JuliaLang/julia · GitHub so that your issue can help the developers as much as possible to reproduce and then fix the bug.