Generation of the matrix corresponding to a map composition fails

[andreasvarga]

I recently implemented two operators Dl and De, known in the literature as elimination operator and duplication operator, respectively. The matrices of these operators for 2-by-2 matrices can be constructed as follows:

julia> using LinearAlgebra, LinearMaps
julia> import MatrixEquations: triu2vec, vec2triu
julia> De = matelimop(2);  # use code below
julia> Dl = matduplop(2);  # use code below

julia> Matrix(Dl)
4×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

julia> Matrix(De)
3×4 Matrix{Float64}:
 1.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0
 0.0  0.0  0.0  1.0

I also explicitly implemented the adjunct operators for which, as expected, I obtained

julia> Matrix(Dl')
3×4 Matrix{Float64}:
 1.0  0.0  0.0  0.0
 0.0  1.0  1.0  0.0
 0.0  0.0  0.0  1.0

julia> Matrix(De')
4×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

For the composed map Dl*De I obtain

julia> Matrix(De*Dl)
3×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

but for its adjoint, explicitly computed, the generation of the matrix fails, with a strange error message:

julia> Matrix(Dl'*De')
ERROR: transpose not implemented for 4×3 MatrixEquations.matduplop{Float64}
Stacktrace:
  [1] error(s::String)
    @ Base .\error.jl:35
  [2] _unsafe_mul!(y::SubArray{Float64, 1, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true}, A::LinearMaps.TransposeMap{Float64, MatrixEquations.matduplop{Float64}}, x::Vector{Float64})
    @ LinearMaps C:\Users\Andreas\.julia\packages\LinearMaps\ojziG\src\transpose.jl:54
  [3] _compositemul!(y::SubArray{Float64, 1, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true}, A::LinearMaps.CompositeMap{Float64, Tuple{LinearMaps.TransposeMap{Float64, MatrixEquations.matelimop{Float64}}, LinearMaps.TransposeMap{Float64, MatrixEquations.matduplop{Float64}}}}, x::Vector{Float64}, source::Nothing, dest::Nothing)
    @ LinearMaps C:\Users\Andreas\.julia\packages\LinearMaps\ojziG\src\composition.jl:193
  [4] _compositemul!(y::SubArray{Float64, 1, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true}, A::LinearMaps.CompositeMap{Float64, Tuple{LinearMaps.TransposeMap{Float64, MatrixEquations.matelimop{Float64}}, LinearMaps.TransposeMap{Float64, MatrixEquations.matduplop{Float64}}}}, x::Vector{Float64})
    @ LinearMaps C:\Users\Andreas\.julia\packages\LinearMaps\ojziG\src\composition.jl:191
  [5] _unsafe_mul!(y::SubArray{Float64, 1, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true}, A::LinearMaps.CompositeMap{Float64, Tuple{LinearMaps.TransposeMap{Float64, MatrixEquations.matelimop{Float64}}, LinearMaps.TransposeMap{Float64, MatrixEquations.matduplop{Float64}}}}, x::Vector{Float64})
    @ LinearMaps C:\Users\Andreas\.julia\packages\LinearMaps\ojziG\src\composition.jl:177
  [6] mul!
    @ C:\Users\Andreas\.julia\packages\LinearMaps\ojziG\src\LinearMaps.jl:173 [inlined]
  [7] _generic_map_mul!(Y::Matrix{Float64}, A::LinearMaps.CompositeMap{Float64, Tuple{LinearMaps.TransposeMap{Float64, MatrixEquations.matelimop{Float64}}, LinearMaps.TransposeMap{Float64, MatrixEquations.matduplop{Float64}}}}, s::Bool)
    @ LinearMaps C:\Users\Andreas\.julia\packages\LinearMaps\ojziG\src\LinearMaps.jl:323
  [8] _unsafe_mul!
    @ C:\Users\Andreas\.julia\packages\LinearMaps\ojziG\src\LinearMaps.jl:275 [inlined]
  [9] Matrix{Float64}(A::LinearMaps.CompositeMap{Float64, Tuple{LinearMaps.TransposeMap{Float64, MatrixEquations.matelimop{Float64}}, LinearMaps.TransposeMap{Float64, MatrixEquations.matduplop{Float64}}}})
    @ LinearMaps C:\Users\Andreas\.julia\packages\LinearMaps\ojziG\src\conversion.jl:4
 [10] (Matrix)(A::LinearMaps.CompositeMap{Float64, Tuple{LinearMaps.TransposeMap{Float64, MatrixEquations.matelimop{Float64}}, LinearMaps.TransposeMap{Float64, MatrixEquations.matduplop{Float64}}}})
    @ LinearMaps C:\Users\Andreas\.julia\packages\LinearMaps\ojziG\src\conversion.jl:12
 [11] top-level scope
    @ REPL[76]:1

I wonder if the error is mine or perhaps is a problem of LinearMaps.

The code to define the above operators is given below:

"""
    M = matelimop(n)

Define the elimination operator of all `n×n` matrices to select their upper triangular parts.  
See [here](https://en.wikipedia.org/wiki/Duplication_and_elimination_matrices) for the 
definition of an elimination matrix for the selection of lower triangular parts.
"""
struct matelimop{T} <: LinearMaps.LinearMap{T}
   dim::Int
   function matelimop{T}(dim::Int) where {T}
      dim >= 0 || throw(ArgumentError("dimension must be non-negative"))
      return new{T}(dim)
   end
   function matelimop(n::Int)
      (n ≥ (0))  || throw(ArgumentError("dimension must be non-negative"))
      return new{Int}(n)
   end
end
Base.size(A::matelimop) = (n = A.dim; return (div(n*(n+1),2), n*n))

function mul!(y::AbstractVector, A::matelimop, x::AbstractVector)
   X = reshape(x, A.dim, A.dim)
   LinearMaps.check_dim_mul(y, A, x)
   copyto!(y, triu2vec(X))
   return y
end
function mul!(x::AbstractVector, AT::LinearMaps.TransposeMap{T,<:matelimop{T}}, y::AbstractVector) where {T}
   Y = vec2triu(y)
   LinearMaps.check_dim_mul(x, AT, y)
   copyto!(x, Y[:])
   return x
end

matelimop(A) = matelimop{eltype(A)}(size(A,1))

"""
    M = matduplop(n)

Define the duplication operator of all `n×n` matrices to reconstruct a hermitian matrix from its upper triangular elements.
See [here](https://en.wikipedia.org/wiki/Duplication_and_elimination_matrices) for the 
definition of a duplication matrix from the lower triangular parts.
"""
struct matduplop{T} <: LinearMaps.LinearMap{T}
   dim::Int
   function matduplop{T}(dim::Int) where {T}
      dim >= 0 || throw(ArgumentError("dimension must be non-negative"))
      return new{T}(dim)
   end
   function matduplop(n::Int)
      (n ≥ (0))  || throw(ArgumentError("dimension must be non-negative"))
      return new{Int}(n)
   end
end
Base.size(A::matduplop) = (n = A.dim; return (n*n, div(n*(n+1),2)))

function mul!(y::AbstractVector, A::matduplop, x::AbstractVector)
   y[:] = vec2triu(x, her = true)
   # LinearMaps.check_dim_mul(y, A, x)
   # copyto!(y, triu2vec(X))
   return y
end
function mul!(x::AbstractVector, AT::LinearMaps.TransposeMap{T,<:matduplop{T}}, y::AbstractVector) where {T}
   Y = reshape(y, (AT.lmap.dim, AT.lmap.dim))
   LinearMaps.check_dim_mul(x, AT, y)
   copyto!(x, triu2vec(Y+Y'-Diagonal(Y)))
   return x
end

matduplop(A) = matduplop{eltype(A)}(size(A,1))

[andreasvarga]

I fixed this issue by replacing mul! with LinearMaps._unsafe_mul! in the above definitions of methods for transposes. I wonder if this is the only way to fix this issue.

[dkarrasch]

Sorry for my lack of response. That might be the right fix. For most cases, there's a fallback for _unsafe_mul! that calls mul!, but when the map is wrapped in a TransposeMap, that doesn't seem to work. Perhaps we should state more clearly that it is actually _unsafe_mul! that should be overloaded. In reward, one doesn't have to do size checking yourself.

[dkarrasch]

I realized that we actually do recommend to implement LinearMaps._unsafe_mul! for custom map types in our documentation: https://julialinearalgebra.github.io/LinearMaps.jl/stable/generated/custom/. Given those, the usual LinearAlgebra.mul! interface just works.

As for the action of those operators, can't you make them in-place functions? Like

_unsafe_mul!(y::AbstractVector, A::matelimop, x::AbstractVector) =
    triu2vec!(y, x)

etc.? You would then avoid intermediate allocations and extra copyto!s.

[andreasvarga]

The actual implementation (with changed operator name) is:

function LinearMaps._unsafe_mul!(y::AbstractVector, A::duplicationop, x::AbstractVector)
   # y[:] = vec2triu(x, her = true)
   n = A.dim
   Y = reshape(y,n,n)
   @inbounds begin
      k = 1
      for j = 1:n
          for i = 1:j
              Y[i,j] = x[k]
              k += 1
          end
          for i = j+1:n
              Y[i,j] = conj(Y[j,i])
          end
      end
   end
end

I think this uses the least memory allocation.

[dkarrasch]

In your case, the reshape is relevant only for the index reassignment. You could avoid it by doing the index manipulation manually:

function LinearMaps._unsafe_mul!(y::AbstractVector, A::duplicationop, x::AbstractVector)
   # y[:] = vec2triu(x, her = true)
   n = A.dim
   @inbounds begin
      k = 1
      for j = 1:n
          for i = 1:j
              y[(i-1)*n + j] = x[k]
              k += 1
          end
          for i = j+1:n
              y[(i-1)*n + j] = conj(y[(j-1)*n + i])
          end
      end
    end
    return y
end

Perhaps it's worth benchmarking.

[andreasvarga]

julia> n = 100;
julia> Dl = duplicationop(n);
julia> x = rand(ComplexF64,Int((n*(n+1)/2))); y = zeros(ComplexF64,n*n);
julia> @time LinearMaps._unsafe_mul!(y, Dl, x);
  0.000022 seconds

Perfect!

[andreasvarga]

I did similarly for the adjoint/transpose:

function LinearMaps._unsafe_mul!(x::AbstractVector, AT::LinearMaps.TransposeMap{T,<:duplicationop{T}}, y::AbstractVector) where {T}
   n = AT.lmap.dim
   #Y = reshape(y, n, n)
   LinearMaps.check_dim_mul(x, AT, y)
   # x[:] = triu2vec(Y+transpose(Y)-Diagonal(Y))
   @inbounds begin
      k = 1
      for j = 1:n
         for i = 1:j
            #x[k] = i == j ? Y[j,j] : Y[i,j] + Y[j,i]
            x[k] = i == j ? y[(j-1)*n + j] : y[(i-1)*n + j] + y[(j-1)*n + i]
            k += 1
         end
      end
   end
end

julia> n = 100;

julia> x = rand(ComplexF64,Int((n*(n+1)/2))); y = zeros(ComplexF64,n*n);

julia> Dl = duplicationop(n);

julia> @time LinearMaps._unsafe_mul!(x, transpose(Dl), y);
  0.000024 seconds (1 allocation: 16 bytes)

[dkarrasch]

You can also skip the check_dim_mul. That is already taken care of by the generic mul! implementation in LinearMaps.jl. Just have the naked code (+ potential require_one_based_indexing if necessary).

[andreasvarga]

Thanks for the hint. I implemented similarly the operations for the elimination operator.

Actually, I don't use these operators in MatrixEquations, but they are very useful for testing purposes. Using them, I was able to fix the transpose/adjoint issues for half-vectorization based Lyapunov operators for real data. The complex case is still not fixed and I wonder if operations like complex conjugation are linear (I guess not), which could be the cause of my problems. I needed to fix the transpose issues, because I started to add several iterative solvers based on the conjugate gradient method and the operations with adjoint/transpose were not correctly implemented for half-vectorization based approaches. At a certain moment. before issuing a new minor release, I would be very pleased if you could have a look on these new developments.

One pleasant issue I would like to mention is that the defined operators works also for sparse matrices and can be also used with the two alternative solvers lsqr and lsmr available in the IterativeSolvers package.

I am a little bit excited using the iterative methods, being able to solve now arbitrary linear matrix equations of the very general form ∑ A_i*X*B_i + ∑ C_j*transpose(X)*D_j = E, where any of matrices can be equal toI. I wonder if it would be possible to extend the definitions of some operators to also cover the case. when the above A_i, B_i, ... are linear maps.

If you consider appropriate, you can close this issue.

[dkarrasch]

I wonder if it would be possible to extend the definitions of some operators to also cover the case. when the above A_i, B_i, ... are linear maps.

It depends on which operations you use to create new operators. If those operations produce a new LinearMap, then that should just work with IterativeSolvers.jl, unless they require AbstractMatrix.