Solving Lyapunov equation for matrices given in Adjoint form

[guberger]

There seems to be no consensus about what should be the Lyapunov equation: A'P+PA+Q=0 or AP+PA'+Q=0. For instance, Wikipedia says differently than the doc of lyap. For this reason, I think it would be nice to have an easy interface to switch between both definitions by allowing to put A' directly in the arguments of lyap.

[simonbyrne]

Is it possible to compute it without first explicitly computing the adjoint?

[simonbyrne]

In any case, this should have a test.

[imciner2]

We always need to form the full matrices to pass into A and C for the lyap function to work, since it is only defined for StridedArray. This behavior seems similar to how other methods such as schur implement their dispatching on the more generic arrays.

[imciner2]

Why not generalize this to be an AbstractMatrix for both A and C:

lyap(A::AbstractMatrix{T}, C::AbstractMatrix{T}) where T = lyap(Matrix(A), Matrix(C))

That will allow for the adjoints to be passed in, as well as other matrix types to also be used (such as Tridiagonal, etc.) and make this bare-bones example code work (it will not currently work because it doesn't allow the Tridiagonal dispatch on the second argument).

A = randn(5,5)
B = randn(5,5)

lyap( A, Tridiagonal( B'*B ) )

[dkarrasch]

This should be

function lyap(A::AbstractMatrix, C::AbstractMatrix)
    T = promote_type(float(eltype(A)), float(eltype(C)))
    return lyap(copy_to_array(A, T), copy_to_array(C, T))
end

This will make sure that the matrices are strided and have float eltype, in one go.