Doc Instructions for Transpose Insufficient for Full Functionality

[miakramer]

I've run into an issue where the instructions on creating transposes for custom linear maps are insufficient for the custom map to fully work. A MWE using the documented example of MyFillMap:

using LinearMaps, LinearAlgebra

struct MyFillMap{T} <: LinearMaps.LinearMap{T}
    λ::T
    size::Dims{2}
    function MyFillMap(λ::T, dims::Dims{2}) where {T}
        all(≥(0), dims) || throw(ArgumentError("dims of MyFillMap must be non-negative"))
        promote_type(T, typeof(λ)) == T || throw(InexactError())
        return new{T}(λ, dims)
    end
end

Base.size(A::MyFillMap) = A.size

function LinearAlgebra.mul!(y::AbstractVecOrMat, A::MyFillMap, x::AbstractVector)
    LinearMaps.check_dim_mul(y, A, x)
    return fill!(y, iszero(A.λ) ? zero(eltype(y)) : A.λ*sum(x))
end

function LinearAlgebra.mul!(
    y::AbstractVecOrMat,
    transA::LinearMaps.TransposeMap{<:Any,<:MyFillMap},
    x::AbstractVector
)
    LinearMaps.check_dim_mul(y, transA, x)
    λ = transA.lmap.λ
    return fill!(y, iszero(λ) ? zero(eltype(y)) : transpose(λ)*sum(x))
end

Then,

julia> let M = MyFillMap(1, (5, 5))' * LinearMaps.UniformScalingMap(1, (5, 5))
         M * zeros(5)
       end
ERROR: transpose not implemented for 5×5 MyFillMap{Int64}
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:33
 [2] _unsafe_mul!(y::Vector{Float64}, A::LinearMaps.TransposeMap{Int64, MyFillMap{Int64}}, x::Vector{Float64})
   @ LinearMaps ~/.julia/packages/LinearMaps/NF1pj/src/transpose.jl:54
 [3] _compositemul!(y::Vector{Float64}, A::LinearMaps.CompositeMap{Int64, Tuple{LinearMaps.UniformScalingMap{Int64}, LinearMaps.TransposeMap{Int64, MyFillMap{Int64}}}}, x::Vector{Float64}, source::Vector{Float64}, dest::Nothing)
   @ LinearMaps ~/.julia/packages/LinearMaps/NF1pj/src/composition.jl:186
 [4] _compositemul!(y::Vector{Float64}, A::LinearMaps.CompositeMap{Int64, Tuple{LinearMaps.UniformScalingMap{Int64}, LinearMaps.TransposeMap{Int64, MyFillMap{Int64}}}}, x::Vector{Float64})
   @ LinearMaps ~/.julia/packages/LinearMaps/NF1pj/src/composition.jl:185
 [5] _unsafe_mul!(y::Vector{Float64}, A::LinearMaps.CompositeMap{Int64, Tuple{LinearMaps.UniformScalingMap{Int64}, LinearMaps.TransposeMap{Int64, MyFillMap{Int64}}}}, x::Vector{Float64})
   @ LinearMaps ~/.julia/packages/LinearMaps/NF1pj/src/composition.jl:168
 [6] mul!(y::Vector{Float64}, A::LinearMaps.CompositeMap{Int64, Tuple{LinearMaps.UniformScalingMap{Int64}, LinearMaps.TransposeMap{Int64, MyFillMap{Int64}}}}, x::Vector{Float64})
   @ LinearMaps ~/.julia/packages/LinearMaps/NF1pj/src/LinearMaps.jl:163
 [7] *(A::LinearMaps.CompositeMap{Int64, Tuple{LinearMaps.UniformScalingMap{Int64}, LinearMaps.TransposeMap{Int64, MyFillMap{Int64}}}}, x::Vector{Float64})
   @ LinearMaps ~/.julia/packages/LinearMaps/NF1pj/src/LinearMaps.jl:131
 [8] top-level scope
   @ REPL[5]:2

Adding this:

function LinearMaps._unsafe_mul!(
    y::AbstractVecOrMat,
    transA::LinearMaps.TransposeMap{<:Any,<:MyFillMap},
    x::AbstractVector
)
    λ = transA.lmap.λ
    fill!(y, iszero(λ) ? zero(eltype(y)) : transpose(λ)*sum(x))
end

Makes the example work.

[dkarrasch]

Thanks for your report. I noticed this issue in a different context, therefore I left a comment

LinearMaps.jl/src/LinearMaps.jl

Lines 201 to 202 in 98843dd

	function _generic_mapvec_mul!(y, A, x, α, β)
	# this function needs to call mul! for, e.g., AdjointMap{...,<:CustomMap}

and a couple of tests

LinearMaps.jl/test/linearmaps.jl

Line 106 in 98843dd

@testset "transpose of new LinearMap type" begin

but only recently realized that this will still fail in "higher-order" linear maps. I believe the additions in
https://github.com/JuliaLinearAlgebra/LinearMaps.jl/pull/173/files#diff-c8422fd2fc31c8d6e3ba2804f0ba3aa0a152cb50a8d7481514e891abb82dfa94

will solve the issue here. In particular, we will start to recommend to always write LinearMaps._unsafe_mul! methods.