search

JuliaLinearAlgebra / LinearMaps.jl

A Julia package for defining and working with linear maps, also known as linear transformations or linear operators acting on vectors. The only requirement for a LinearMap is that it can act on a vector (by multiplication) efficiently.
Other
303 stars 42 forks source link

Doc Instructions for Transpose Insufficient for Full Functionality #181

Closed miakramer closed 2 years ago

miakramer commented 2 years ago

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:

```julia 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 commented 2 years ago

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

https://github.com/JuliaLinearAlgebra/LinearMaps.jl/blob/98843dd43adc5788045d6f8ccdc703a8a33ae653/src/LinearMaps.jl#L201-L202

and a couple of tests

https://github.com/JuliaLinearAlgebra/LinearMaps.jl/blob/98843dd43adc5788045d6f8ccdc703a8a33ae653/test/linearmaps.jl#L106

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.