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)

[dkarrasch]

I've been thinking and working on that issue here now for a while. It seems to me, that this is probably not an issue with LinearMaps.jl. Both attempts to solve it lead to value-dependent branches of code, or even to type instability altogether. But Julia has a type systems, and dispatch should (in a type-stable environment) should do the job.

Given that A * I is not a no-op for A::AbstractMatrix, I don't quite understand why making ::LinearMap + I return the linear map would solve the problem with superfluous operations. You would still have extra effort in the matrix case. So, IMHO, I think this issue should be handled on the downstream package side. For instance, if there already is an Identity type, then make generically

(*)(A, ::Identity) = A
(*)(::Identity, A) = A

Now, in some linalg function, you could have

function myfun(A, args...; P=Identity(), kwargs...)
    some_stuff(A*P)
    ....
end

where you could then even do dispatch on the type of A and what not. So, for the user, using "no preconditioner" would not mean to supply I, but to leave empty the kwarg P. Well, or have P as a positional arg with default value Identity. At least, I don't seem to find any indication in the ecosystem that, in multiplication with a UniformScaling, anyone would check the value of λ and change behavior based on that.

Generally, I came to appreciate the purity of how this package was designed by @Jutho. As an example, you can have a real matrix, but wrap it in a LinearMap{ComplexF64} type. That doesn't mean that the package is going to promote the underlying array to that complex type. The same is with, say, addition of uniform scalings. The resulting LinearCombination gets the promoted eltype, but the λ itself of the uniform scaling is not touched. In this spirit, I believe if someone provides a uniform scaling where λ=1 or λ = true happens to be, then this should be respected and performed (as is done in Base with matrices). If the user wants something to act like nothing in multiplication, I think this calls for a type that is defined as such, and not for a re-interpretation of a scaling object.

Sorry for the long text, but I suffered quite some pain while trying and failing, just to realize that I finally agree with my initial gut feeling.

PS: Looking at some of the earlier comments, there was concern that we "create y's" just to copy to. Don't worry! We don't create unnecessary ys. What I was describing is how multiplication with a unit uniform scaling (not necessarily of boolean type) works in mul!(y, ::UniformScaling, x): if λ==1, we don't even multiply x with 1, but outright copy x to y, which also handles possible type promotion (again, without performing the element-wise multiplication). As you can quickly test, this is much faster than doing the trivial multiplication.