For a custom map L Matrix(L-L) is not a zero matrix

[andreasvarga]

I implemented a custom linear map for the Lyapunov operator L(X) = A*X+X*A', where A and X are nxn matrices, In the simple case when n = 1and A = 1, this is just the linear function L(X) = 2*X. Here is how the Lyapunov operator L is defined for this simple case and then applied to a vector with a single component:

julia> L=lyapop(1)
1×1 lyapop{Int64}

julia> L*[1]
1-element Vector{Int64}:
 2

julia> -L*[1]
1-element Vector{Int64}:
 -2

Also, the corresponding matrix representation is correctly determined:

julia> Matrix(L)
1×1 Matrix{Int64}:
 2

The matrix representation of the operator L-L should be zero, but I am getting the wrong result

julia> Matrix(L-L)
1×1 Matrix{Int64}:
 4

and further some other nonsense results:

julia> Matrix(-L)
1×1 Matrix{Int64}:
 2

julia> Matrix(0.1*L)
1×1 Matrix{Float64}:
 2.0

julia> Matrix(10*L)
1×1 Matrix{Int64}:
 2

so the scaling ofL is completely ignored.

I defined a single mul! function of the form mul!(y, L, x) , which apparently does it job, but I am not sure if I have to look for an error in my code or is an error somewhere at the level of composing operators, since

julia> -L
1×1 LinearMaps.ScaledMap{Int64} with scale: -1 of
  1×1 lyapop{Int64}

I would appreciate any hint to solve this problem and I can, of course, provide further details of my implementation.

[dkarrasch]

Thanks for the report. I think I'd need to see the definition of the type and the mul! method. BTW, why do you need your own custom type. I'd define the Lyaponov operator via the Kronecker sum, kronsum(A, A') or simply A ⊕ A'.

julia> using LinearMaps

julia> A = fill(1, 1, 1)
1×1 Matrix{Int64}:
 1

julia> L = kronsum(A, A')
1×1 LinearMaps.KroneckerSumMap{Int64} with 2 maps:
  1×1 LinearMaps.WrappedMap{Int64} of
    1×1 Matrix{Int64}
  1×1 LinearMaps.WrappedMap{Int64} of
    1×1 adjoint(::Matrix{Int64}) with eltype Int64

julia> Matrix(L)
1×1 Matrix{Int64}:
 2

julia> Matrix(-L)
1×1 Matrix{Int64}:
 -2

julia> Matrix(10L)
1×1 Matrix{Int64}:
 20

julia> Matrix(0.1L)
1×1 Matrix{Float64}:
 0.2

[andreasvarga]

Here is the definition of the Lyapunov operator, where I commented out the check for Schur form:

"""
    L = lyapop(A; disc = false, her = false) 

Define, for an `n x n` matrix `A`, the continuous Lyapunov operator `L:X -> AX+XA'`
if `disc = false` or the discrete Lyapunov operator `L:X -> AXA'-X` if `disc = true`.
If `her = false` the Lyapunov operator `L:X -> Y` maps general square matrices `X`
into general square matrices `Y`, and the associated matrix `M = Matrix(L)` is 
``n^2 \\times n^2``.
If `her = true` the Lyapunov operator `L:X -> Y` maps symmetric/Hermitian matrices `X`
into symmetric/Hermitian matrices `Y`, and the associated matrix `M = Matrix(L)` is 
``n(n+1)/2 \\times n(n+1)/2``.
For the definitions of the Lyapunov operators see:

M. Konstantinov, V. Mehrmann, P. Petkov. On properties of Sylvester and Lyapunov
operators. Linear Algebra and its Applications 312:35–71, 2000.
"""
struct lyapop{T} <: LinearMaps.LinearMap{T}
  A::AbstractMatrix
  disc::Bool
  her::Bool
  adj::Bool
  sf::Bool
  function lyapop(A::AbstractMatrix{T}; disc = false, her = false) where {T}
      LinearAlgebra.checksquare(A)
      adj = isa(A,Adjoint)
      # schur_flag = adj ? isschur(A.parent) : isschur(A)
      schur_flag = false
      return new{T}(adj ? A.parent : A, disc, her, adj, schur_flag)
  end
end
lyapop(A::Number; kwargs...) = lyapop(A*ones(eltype(A),1,1); kwargs...) 
function Base.size(L::lyapop)
    n = size(L.A,1)
    N = L.her ? Int(n*(n+1)/2) : n*n
    return (N,N)
end
LinearAlgebra.adjoint(L::lyapop{T}) where T = lyapop(L.adj ? L.A : L.A'; disc = L.disc, her = L.her)
LinearAlgebra.transpose(L::lyapop{T}) where T = lyapop(L.adj ? L.A : L.A'; disc = L.disc, her = L.her)

And here is themul! function in a simpler form to allow execution for you:

function LinearAlgebra.mul!(y::AbstractVector, L::lyapop{T}, x::AbstractVector) where T
  n = size(L.A,1)
  T1 = promote_type(T, eltype(x))
  if L.her
    #X = vec2triu(convert(Vector{T1}, x), her = true)
   # if L.disc
   #   L.adj ? (y[:] = triu2vec(utqu(X,L.A) - X)) : (y[:] = triu2vec(utqu(X,L.A') - X))
   # else
   #    L.adj ? Y = L.A' * X : Y = L.A * X
   #    y[:] = triu2vec(Y + Y')
   #end
  else
    X = reshape(convert(Vector{T1}, x), n, n)
    if L.disc
       L.adj ? (y[:] = (L.A'*X*L.A - X)[:]) : (y[:] = (L.A*X*L.A' - X)[:])
    else
       L.adj ? (y[:] = (L.A'*X + X*L.A)[:]) : (y[:] = (L.A*X + X*L.A')[:])
    end
  end
  return y[:]
end

As you can see, the definition is more involved and the Kronecker-based definition does not cover the symmetric mapping case. Moreover, the definition of the inverse Lyapunov mapping involves the use of special solvers for Lyapunov and Sylvester matrix equations (e.g., as those available in MatrixEquations.jl).

[andreasvarga]

I enhanced the mul! function to reduce as much as possible the memory allocations. The execution times are comparable with the Kronecker-sum based computation, which is to be expected. However, the memory allocation is the same as with the Kronecker-sum, for which I have no explanation. I can just speculate, but in my opinion internally I am storing in lyapop a 500x500 matrix, which accounts for 1.907 MiB memory, while building the Kronecker-sum should involve the double of this (i.e., 3.82 MiB). These are the benchmark results I obtained:

julia> A = rand(500,500);

julia> x = rand(250000);

julia> LK = kronsum(A,A);

julia> L = lyapop(A);

julia> norm((L-L)*x)   # this line I entered erroneously (by omitting K in LK, see below)
0.0

julia> norm((L-LK)*x)
0.0

julia> @btime LK*x;
  4.748 ms (7 allocations: 1.91 MiB)

julia> @btime L*x;
  4.925 ms (17 allocations: 1.91 MiB)

julia> @btime LK'*x;
  4.924 ms (8 allocations: 1.91 MiB)

julia> @btime L'*x;
  4.951 ms (20 allocations: 1.91 MiB)

For your information, I replaced in mul!

L.adj ? (y[:] = (L.A'*X + X*L.A)[:]) : (y[:] = (L.A*X + X*L.A')[:])

with

        Y = reshape(y, (n, n))   
        if L.adj 
           #(y[:] = (L.A'*X + X*L.A)[:]) 
           mul!(Y, L.A', X)
           mul!(Y, X, L.A, true, true)
        else
           #(y[:] = (L.A*X + X*L.A')[:])
           mul!(Y, X, L.A')
           mul!(Y, L.A, X, true, true)
        end

and removedreturn y[:]. The error is gone!

julia> L = lyapop(1);

julia> Matrix(L)
1×1 Matrix{Int64}:
 2

julia> Matrix(-L)
1×1 Matrix{Int64}:
 -2

julia> Matrix(10L)
1×1 Matrix{Int64}:
 20

julia> Matrix(0.1L)
1×1 Matrix{Float64}:
 0.2

Perhaps you have an explanation for this. Many thanks for your time.

[dkarrasch]

Awesome! I was going to suggest that you copy the kronsum code to avoid the allocations in the one branch. I don't use the y[:] pattern at all, so I would have needed to go and see what exactly it does, like allocating a new object, shared memory etc.

I'm not sure I understand why you want to store the matrix by all means, i.e., stripping off the Adjoint wrapper, when you then have two formulae that do nothing but to account for a possible Adjoint. Unless you have other reasons, I think you could safely store any matrix and apply one formula.

[dkarrasch]

BTW, it should be LK = kronsum(A',A);. It might change runtimes a bit due to memory-layout.

[dkarrasch]

Also, the allocations are almost entirely due to the final result:

julia> A = rand(500,500);

julia> LK = kronsum(A',A);

julia> x = rand(250000);

julia> y = similar(x);

julia> @btime mul!($y, $LK, $x);
  3.512 ms (5 allocations: 208 bytes)

julia> @btime $LK * $x;
  3.699 ms (7 allocations: 1.91 MiB)

The few bytes left are for the "reshape"-wrapper, and probably unavoidable, or at least harmless.

[andreasvarga]

I think LK = kronsum(A,A) is correct, because(A ⊗ I + I ⊗ A)vec(X) corresponds to A*X+X*A' and kronsum(A',A') corresponds to the adjoint operator A'*X+X*A.

Storing one matrix for both A and A' is relevant, when A is in a Schur form (you are right, this plays no role in the case for lyapop). But for the inverse mapping, this information is explicitly exploited by the solvers, which, depending on the information on adjoint, use one algorithm or another. In this way, all algorithms in MatrixEquations work with the upper Schur form, as it results from standard computations with the function schur.

I was not aware that usingreturn y[:] adds the whole memory occupied by y. Could this be also the cause that -L was not correctly handled ? I will check all my software to see if this issue caused unnecessary memory allocations.

Many thanks for this stimulating discussion. If you consider appropriate, you can close this issue.

[dkarrasch]

Are your matrices always real? You are right about kronsum(A, A), because by the vec-trick, (A ⊗ I + I ⊗ A)vec(X) corresponds to A*X+X*transpose(A). So in the complex case, there may be issues?

[andreasvarga]

Yes, it should be kronsum(transpose(A),transpose(A)) in the second case. And yes, I cover all BlasFloat cases.

[andreasvarga]

I was wrong in my previous comment. The complex case needs a special handling when using Kronecker products.

Daniel Karrasch @.***> schrieb am Sa., 26. Juni 2021, 17:13:

Are your matrices always real? You are right about kronsum(A, A), because by the vec-trick, (A ⊗ I + I ⊗ A)vec(X) corresponds to AX+Xtranspose(A). So in the complex case, there may be issues?

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/Jutho/LinearMaps.jl/issues/155#issuecomment-869016795, or unsubscribe https://github.com/notifications/unsubscribe-auth/ALJDHED4D6QJJDV2AOSOLYDTUXVDLANCNFSM47KALMPA .

[dkarrasch]

I guess it should be something like kronsum(A, conj(A)), right? Which you could use as a general rule, because conj(A) is a noop for real matrices.

I'm glad you solved the problem and you're using our package. I'll close the issue, but you're most welcome to report any other issues you encounter, of course.

[andreasvarga]

Just a final remark on why to prefer explicit definition like L = oplyap(A) for Lyapunov mapping instead of LK = kronsum(A,A).

If you dig a lower in the resulting objects L and LK, you see that L stores just a single copy of A, while the LK object stores A twice, which certainly is not desirable in this case:

julia> A = rand(500,500);

julia> L = oplyap(A);  # in the meantime I changed the name from lyapop to oplyap! 

julia> LK = kronsum(A,A);

julia> fieldnames(typeof(L))
(:A, :disc, :her, :adj, :sf)

julia> size(L.A)
(500, 500)

julia> fieldnames(typeof(LK.maps[1]))
(:lmap, :_issymmetric, :_ishermitian, :_isposdef)

julia> size(LK.maps[1].lmap)
(500, 500)

julia> size(LK.maps[2].lmap)
(500, 500)

In my opinion, the memory allocation issue in LM is somehow hidden. A function which counts all underlying data to a map defintion would be useful to globally analyse the memory usage, since this information is not provided by @time:

julia> @time L = oplyap(A);
  0.000005 seconds (1 allocation: 32 bytes)

julia> @time LK = kronsum(A,A);
  0.000012 seconds (1 allocation: 48 bytes)

[dkarrasch]

No, there is no "second copy" of the matrix. It's true that kronsum stores two copies of the reference to A, no matter if A is a matrix or another linear map. In memory, there exists only one copy of the matrix. And this is what you see in the @time macro: the extra bytes are there to store the second reference.

[andreasvarga]

Thanks. These are good news. Then, I wonder if it is possible to rely on kronsum in the defintion of oplyap, in the case when the hermitian mapping option is not used. Should I possibly need to define two maps, i.e., one for each case?

[dkarrasch]

Yes, instead of copying the mul! code, you could construct a KronSumMap in your mul! method. Roughly like this:

LinearAlgebra.mul!(y::AbstractVector, A::oplyap, x::AbstractVector)
    if herm
        ...
    else
        mul!(y, kronsum(A.A, A.A), x)
    end
    return y
end

The kronsum call itself only calls the lazy constructor, so this should take less than a microsecond:

julia> using LinearMaps

julia> using BenchmarkTools

julia> A = rand(500, 500);

julia> fun(A) = kronsum(A, A)
fun (generic function with 1 method)

julia> @btime fun($A);
  524.269 ns (0 allocations: 0 bytes)

It depends on how much you want to avoid a little code duplication and how much overhead you obtain relative to the multiplication cost.

[andreasvarga]

It seem that it is possible to use a single definition. Thanks.

[andreasvarga]

The following (only thedisc = false case implemented) also works:

struct oplyap{T} <: LinearMaps.LinearMap{T}
  A::AbstractMatrix
  disc::Bool
  her::Bool
  adj::Bool
  sf::Bool
  function oplyap(A::AbstractMatrix{T}; disc = false, her = false) where {T}
      LinearAlgebra.checksquare(A)
      her || (return kronsum(A,conj(A)))
      adj = isa(A,Adjoint)
      schur_flag = adj ? isschur(A.parent) : isschur(A)
      return new{T}(adj ? A.parent : A, disc, her, adj, schur_flag)
  end
end

The resulting operator is either of oplyap type or of KroneckerSumMap type. In evaluating L*x, each would use its ownmul!. I think this could work.

[dkarrasch]

The problem is that this makes the oplyap constructor type-unstable, which may be acceptable if the construction and the application of the operator are separated by a function barrier. But this is perhaps not quite "recommended Julia style".