RestrictionMap/SamplingMap?

[dkarrasch]

I was reading carefully the PyLops paper that @goretkin mentioned in #98, and noticed that there might be some interest in operators that restrict/sample vectors, y = x[I], whose transpose/adjoint can be interpreted as an embedding y[I] .= x. Interestingly, in that paper they show that their sampling operator is basically as fast as the corresponding multiplication with the corresponding sparse matrix (the authors even compare to the multiplication with a dense sampling matrix and seem to make the case that avoiding to store the matrix is of great benefit). I was curious if we can both avoid storing the matrix and have faster performance compared to the sparse matrix approach, as it should be! Of course, this has nothing to do with LinearMaps.jl, but all with array indexing handled by Julia Base, but still, and it's possible, of course.

This is all trivial and I don't want to blow up this package unnecessarily. But if there's interest, this could be a cute addition.

This is somewhat related to #65, but cannot replace it, since we wouldn't want to represent an IndexMap as the product E*A*R of a transpose restriction map/embedding map E, a linear map A, and a restriction map R (to avoid the intermediate "copies"), but have the indexing included. Thumbs up or down for adding it?

[dkarrasch]

Actually, this would play quite nicely with IndexMap (pseudocode):

struct RestrictionMap <: LinearMap{Bool}
    ...
end
struct EmbeddingMap <: LinearMap{Bool}
    ...
end

transpose(R::RestrictionMap) = EmbeddingMap
transpose(R::EmbeddingMap) = RestrictionMap

Base.:(*)(E::EmbeddingMap, A::LinearMap, R::RestrictionMap) =
    IndexMap(A, (size(E, 1), size(R, 2)), E.indices, R.indices)
Base.:(*)(A::LinearMap, R::RestrictionMap) =
    IndexMap(A, (size(A, 1), size(R, 2)), 1:size(A, 1), R.indices)
Base.:(*)(E::EmbeddingMap, A::LinearMap) =
    IndexMap(A, (size(E, 1), size(A, 2)), E.indices, 1:size(A, 2))

[dkarrasch]

@JeffFessler Do you think there's value in this in LinearMaps.jl? Not sure which other application cases might use it.

[JeffFessler]

This kind of sampling is exactly what I use for MRI reconstruction:
https://github.com/JeffFessler/MIRT.jl/blob/master/src/system/Afft.jl
It is also useful for image inpainting problems.
I put some effort into trying to make it efficient:
https://discourse.julialang.org/t/efficient-non-allocating-in-place-getindex-for-bitarray/42268
The community feedback led to this version of getindex! that I wish were part of Base:
https://github.com/JeffFessler/MIRT.jl/blob/88d3c82ba6f432fa6833c78ee9b2fb4fbf37889d/src/utility/mask.jl
It would be great to have another pair of eyes on this to see if there are even more efficient ways.
I think the "value added" of having it here in LM would be if there is some synergism (e.g., memory savings) by integrating the indexing with the forward op. Just overloading * would be a "convenience" but maybe of minor value, whereas something with mul! might have more impact. Actually you'll notice that the forw! op for my Afft uses a "work" array, which is related to our previous discussion of preallocating work spaces :)
In short, I personally probably would not use a version that provided only * because I already have a more efficient version, but if LM could support mul! in a general way then definitely I would use it.

[dkarrasch]

Thanks for your feedback. You usage with the work array is exactly as it should be: in a specific package that has well-defined context. In general, we may not know what the type of that array should be, right? You could multiply a real map/matrix with a complex vector, so I doubt we want to get into the trouble of defining a default behavior here.

[dkarrasch]

@JeffFessler I only now fully realized that what you're wishing for is addressed by your IndexMap approach, right? I agree that having new types for sampling/embedding alone is not so helpful, unless it melts together with the actual multiplication code as in IndexMap. So I'll close this one here. Do you mind if I take over your PR? Some parts need rewriting due to the changes to the library in the past 3 years since your PR.

[JeffFessler]

The original motivation for IndexMap was to try to simplify and generalize BlockMap, if I recall, and I never got back to it because I guess I didn't really need the generalization and I kind of got lost in 5-arg mul! space. I wasn't thinking about sampling/restriction and it seems likely that those are more compelling uses.
I would be delighted for you to take over the PR and salvage anything useful there!
Thanks!