Make LinearMap * I returns same LinearMap

[platawiec]

Currently,

*(A::LinearMap, I::UniformScaling{Bool})::CompositeMap

That is, a new map is constructed. Could we simply return A in this case (and related cases involving multiplication/division by the identity?)

[JeffFessler]

A complication is that false*I is also of the same type but is mathematically equivalent to 0I.
So a slightly more complicated function would be needed to check the value first.

[dkarrasch]

@JeffFessler That’s correct, which makes me think that A * (::UniformScaling{bool}) would become type-unstable: depending on whether the lambda is true or false, we'd return an object of type of A or a CompositeMap.

@platawiec What exactly is your concern about the CompositeMap? Multiplying with true*I (which is just I) is as expensive as copy!(y, x), see the corresponding mul-code for UniformScalingMaps, so shouldn’t be too expensive. The slightly more interesting case is actually the one @JeffFessler mentioned: when I is actually the zero mapping, then you may not wish to compute the (potentially expensive) A*zero(...), because you already know that it’s going to be zero. But I’m not sure we will want to provide some automatic „tidying up the linear maps“ for the user, but rather compute what the user wants to compute, be it A*I or A*ZeroMap. But knowing the actual concern (or the use case) may help.

[platawiec]

@dkarrasch The specific use-case is for optionally providing pre-conditioners to matrix-free operator representations of linear maps. For instance, IterativeSolvers.jl handles this by supplying a custom Identity type and then dispatching as special-case to reduce the number of ops. I agree that relative to the whatever downstream process you have the call to copyto! is likely negligible, but in memory-constrained scenarios it could be beneficial, where we want to avoid the creation of a y to copyto! in the specific (and common) case of no pre-conditioner.

I suppose the larger point here is how abstractions (pre-conditioners for this particular use-case) and the interface exposed to the user interact with implementation details. I had the idea that it may be appropriate for LinearMaps to handle the "tidying up," as you called it. But I think it's also perfectly valid for "tidying" to lie outside the scope of LinearMaps. The only remaining concern I have is that, if "tidying" is not the concern of LinearMaps, then there may be repeated code in downstream packages, trying to handle the tidying on their own.

[chriscoey]

I would like to avoid copyto! anywhere possible. Maybe there's a different way to avoid that, like using the generic 5-arg mul! perhaps? #56

[JeffFessler]

I was curious what other objects do with * I so I ran a simple test:
Diagonal(1:5) * I returns type Diagonal{Int64,StepRange{Int64,Int64}}
whereas
Diagonal(1:5) has type Diagonal{Int64,UnitRange{Int64}}
which is almost the same type but not quite I guess.
So it seems that even Diagonal times I does not just return the original object, FWIW.
Even worse, Diagonal(1:5) * (0I) throws an error!

[dkarrasch]

Well, as expected, if I introduce

function Base.:(*)(A1::LinearMap, A2::UniformScaling{Bool})
    if A2.λ
        return A1
    else
        return A1 * A2.λ # this could be further reduced to a new ZeroMap type
    end
end

(and similarly the other one) that product becomes type-unstable. The inferred type, however, is a Union of two concrete types:

julia> @inferred A*I
ERROR: return type LinearMaps.WrappedMap{Float64,Array{Float64,2}} does not match i
nferred return type Union{LinearMaps.CompositeMap{Float64,Tuple{LinearMaps.UniformS
calingMap{Bool},LinearMaps.WrappedMap{Float64,Array{Float64,2}}}}, LinearMaps.Wrapp
edMap{Float64,Array{Float64,2}}}

This may be not too bad, actually; at least, AFAIU, Julia is able to handle type instabilities with small type unions. @Jutho any opinion on that matter?

Anyway, I'm not sure we would want to iterate over longer CompositeMaps to see if somewhere in between, there is a zero map, and then iteratively reduce the whole CompositeMap to a ZeroMap.

[JeffFessler]

what about the more general case of, say, 1.0I

function Base.:(*)(A1::LinearMap, A2::UniformScaling) # needs refinement to exclude Bool?
    if A2.λ == 1 # or true ?
        return A1
    else
        return A1 * A2.λ # this could be further reduced to a new ZeroMap type
    end
end

[dkarrasch]

I've been thinking about it. The issue is that this introduces type instability. If we restrict this to only Bool, then the effect is somewhat limited. OTOH, maybe the product with UniformScalings is a special enough case. For instance, 2*L and L*2 are handled by different, type stable methods. This is just the first time we introduce type instability in the package, hoping for union splitting or function barriers in users' code to help us out, that's why I'd like to first hear @Jutho's opinion. Maybe we should have the discussion of the implementation details (and whether we want to do this in the first place) over at #63? You could leave a quick review comment over there, @JeffFessler.

[tkf]

Just FYI, you can use const Id = Init(*) using Init from https://github.com/tkf/InitialValues.jl instead of I to do Id * A === A and A * Id === A.

[JeffFessler]

@platawiec fyi I have added the "simply return A" feature to https://github.com/JeffFessler/LinearMapsAA.jl
(master branch only, not yet tagged)