Add slicing functionality

[dkarrasch]

This is an alternative approach to #145, without adding yet another LinearMap subtype. As it stands now, there is no switch for turning on/off indexing/slicing. I was hoping one could have a submodule and "hide" it from the user, unless s/he asks for using LinearMaps.GetIndex, but that doesn't seem to work as easy as I was hoping for. OTOH, we might as well trust our users that they are cautious about indexing/slicing and not put additional stones in their way.

This may require a few more tests to keep coverage high, some announcement in the docs and some information on how to extend that functionality to own map types.

Closes #145, closes #38.

[codecov]

Codecov Report

Merging #165 (93270b6) into master (7363170) will increase coverage by 0.02%.
The diff coverage is 100.00%.

@@            Coverage Diff             @@
##           master     #165      +/-   ##
==========================================
+ Coverage   99.63%   99.65%   +0.02%     
==========================================
  Files          14       15       +1     
  Lines        1089     1167      +78     
==========================================
+ Hits         1085     1163      +78     
  Misses          4        4              

Impacted Files	Coverage Δ
src/LinearMaps.jl	100.00% <ø> (ø)
src/getindex.jl	100.00% <100.00%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 7363170...93270b6. Read the comment docs.

[codecov]

Codecov Report

Merging #165 (d378b04) into master (12c13fb) will increase coverage by 0.06%.
The diff coverage is 100.00%.

@@            Coverage Diff             @@
##           master     #165      +/-   ##
==========================================
+ Coverage   98.54%   98.61%   +0.06%     
==========================================
  Files          15       16       +1     
  Lines        1307     1372      +65     
==========================================
+ Hits         1288     1353      +65     
  Misses         19       19              

Impacted Files	Coverage Δ
src/LinearMaps.jl	100.00% <ø> (ø)
src/getindex.jl	100.00% <100.00%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 12c13fb...d378b04. Read the comment docs.

[dkarrasch]

Anyone interested in indexing and slicing LinearMaps? 😛 @oschulz

@JeffFessler: if I forward indexing LinearMapsA* to the wrapped LinearMap, LinearMapsAA.jl's tests pass for me locally. Could you please double-check?

If this looks good, we could get this in v3.7 soon, otherwise I'd go and make a release with the other currently merged PRs.

[oschulz]

Anyone interested in indexing and slicing LinearMaps?

Very much so! With this, could we do LinearMap <: AbstractMatrix (#180) as well, to help with dispatch for "non-LinearMap-aware" code?

[dkarrasch]

I'll comment on that in the other issue. For the time being, I'd be interested in getting feedback on the implementation. ;-)

[oschulz]

Ok ;-)

[JeffFessler]

Wow, this is a very comprehensive implementation with several features I didn't know about like Base.Slice and diagind (that therefore are not supported in LinearMapsAA).

I'm not exactly sure how to test this independently.

When you say "forward indexing" do you mean something like this?
Base.getindex(A::LinearMapAX, args...) = Base.getindex(A._lmap, args...)

[dkarrasch]

When you say "forward indexing" do you mean something like this?
Base.getindex(A::LinearMapAX, args...) = Base.getindex(A._lmap, args...)

Exactly.

[JeffFessler]

Personally, I would prefer that anything throws that's "inefficient by construction", so to speak. My suggestion would be any single-element/scalar getindex, as well as row-access in the no-adjoint case.

I think I agree, as long as "scalar" means A[i,j] or A[k] - those are the cases that need warnings (if any).
I consider A[:,j] and probably even A[:,J] to be fine (i.e., no warning needed). Let's also keep in mind that Matrix(A) is same as A[:,:] and there has never been any warning for that case and I'd prefer to retain that behavior.
The general case of A[I,J] also reverts to scalar so also could warn.

[oschulz]

I think I agree, as long as "scalar" means A[i,j] or A[k]

Yes, that's what I meant.

I consider A[:,j] and probably even A[:,J] to be fine (i.e., no warning needed).

I fully agree.

that Matrix(A) is same as A[:,:] and there has never been any warning

Absolutely - that's an efficient operation, after all.

My inuition would be (with scalar indices i, j and array/range-indices I,J) to

• Allow A[:,j], A[:,J], A[I,J] and A[:,:] in general.
• Only allow A[i,:] for maps that support adjoint (need to be careful with complex numbers)
• Disallow A[i,j] in general.

A[I,J] I would allow in general, I think, to support use cases like A[begin+1:end-1:begin+1:end-1] and similar, which would typically still very efficient.

[dkarrasch]

The latest commit largely implements the suggestions (thanks for the clear roadmap!). What I currently have errors on A[i,j], and suggests to consider using A[:,j][i] if one really needs to. Question: shall we do the analogous thing for A[I,j], i.e., error and suggest A[:,j][I]? It allocates the unit vector, the target vector, and then yet another vector after slicing the target. The horizontal slicing behaviour should perhaps be consistent with regard to throwing or working. I wonder what rationale exactly applies to allowing A[I,J] but not A[i,j]? Is it that when we compute entire columns (or rows, whatever is going to require less map applications) and then slice by some vector, that "computing the entire row or column was more worth it" than just for a single component? For the same reason as with A[I,j], A[begin+1:end-1,begin+1:end-1] may not be more efficient than convert(Matrix, A)[begin+1:end-1,begin+1:end-1], especially after #173.

[oschulz]

The latest commit largely implements the suggestions (thanks for the clear roadmap!).

Thanks a lot for doing all of this!

What I currently have errors on A[i,j], and suggests to consider using A[:,j][i] if one really needs to.

Oh, that's a nice approach I think!

Question: shall we do the analogous thing for A[I,j], i.e., error and suggest A[:,j][I]? ... I wonder what rationale exactly applies to allowing A[I,J] but not A[i,j]

My intuition was that code which iterates over A[I,j] with short ranges of I is comparatively rare, while stuff like A[begin+1:end-1,j] is not uncommon (getting the inner part things that were padded to fold with kernels and so on). The same for A[I,J].

A[begin+1:end-1,begin+1:end-1] may not be more efficient than convert(Matrix, A)[begin+1:end-1,begin+1:end-1]

Yes, I think so too. I was assuming that the user code wouldn't want to have the actual matrix in hand, necessarily. Maybe it would in fact be best if A[I,J] would not evaluate anything be return something like a view? But the implementation of A[begin+1:end-1,begin+1:end-1] * x wouldn't be trivial - one would have to pad x first, to the mul, and then do a getindex on the result.

Maybe it's best if we disallow A[I,j], A[i,J] and A[I,J] for now (error with suggestion of alternative code, you you already have) until we come across valid use cases?

[dkarrasch]

Maybe it's best if we disallow A[I,j], A[i,J] and A[I,J] for now (error with suggestion of alternative code, you you already have) until we come across valid use cases?

Yes, I guess that's much better than first allow it and then handle potential "Why is XYZ so slow?" issues... or even disallow it later. I think the main point that the user needs to understand is that these cases do correspond to first compute entire slices and then slice again. If users will have to write that out they will naturally want to think twice if that's really what they want, or whether they were writing it just out of convenience.

[JeffFessler]

Makes sense.
Then probably A[:,:] also should throw for now, whereas Matrix(A) continues to work silently (and possibly slowly)?

[dkarrasch]

Ha, that will make a fantastic PR: basically add a whole bunch of methods that throw! Progress! Yippieh! 🤣

[oschulz]

Then probably A[:,:] also should throw for now,

I think A[:,:] should be allowed. It's very efficient, and it would be surprising for the user if A[:,i] is allowed but A[:,:] isn't.

[JeffFessler]

It's very efficient

I guess that is in the eye of the beholder. If A is the N-point DFT implemented via an FFT (a typical use case for me), then A*x is $O(N \log N)$ whereas A[:,:] is $O(N^2 \log N)$ which is more expensive than the $O(N^2)$ needed to simply construct the matrix in the first place! So my view is that it should be discouraged a bit (if we are also discouraging A[I,J]), but I don't feel too strongly about it.
Perhaps ideally there would be a user switch (like with GPU scalars) but it's not a high enough priority for me to want to try to figure that out.

[oschulz]

then A*x is P(N log N) whereas A[:,:] is O(N^2 log N)

Yes, that's what I mean - it scales proportionally to the number of outputs that you want to produce, both for A[:,i] and A[:,:]. An overhead of log(N) is typically still considered quite efficient in numerics, after all.

I would put it like this: A[:,i] is one mul for an output vector with N numbers - that's definitely Ok, it's the bare performance of the linear map. And A[:,:] is N muls for an output matrix of N^2 numbers. That's very Ok too, I would say.

[dkarrasch]

For higher-dimensional arrays, I guess it is not completely uncommon to put colons in two dimensions. For "matrices", who would construct a (deep) copy via A[:,:], when you have Matrix(A)?

BTW, Matrix(A) was never meant to be "efficient". It was just introduced so that, if you really need (for instance, for stubborn downstream packages 😝 ) to provide a matrix. One often repeated argument by @Jutho is that, for applications that require indexing of many different components, it may be beneficial to pay the full (matrix construction) price once than to have it potentially an unbound number of allocations etc. For such cases, we have Matrix(A) and sparse(A) constructors. So, not "efficient", but perhaps "better than". For that reason, Matrix(A) was never viewed in terms of indexing/slicing. Now, with this PR, these two things (matrix construction and slicing) touch each other via A[:,:], and one may take two perspectives: (a) construction of a matrix, and (b) double-slicing. From the perspective (a), "efficiency" doesn't play a role, it's just supposed to be "better than"; from the perspective (b), this is an awful operation, and should be warned about, to the very least. Obviously, it's just a matter of taste, but I tend to prefer path (b). Note that generic code could use convert(AbstractMatrix, A), which is a noop for AbstractMatrixes.

[oschulz]

For "matrices", who would construct a (deep) copy via A[:,:], when you have Matrix(A)

It's not quite the same in general, because for normal arrays A[:,:] would typically preserve the array type. For LinearMaps it won't matter much which is used, I guess.

My argument for allowing A[:,:] is more that we can easily do it and that allowing A[:,i] (which we want) makes forbidding A[:,:] (which has the same relative efficiency) surprising for the user.

[dkarrasch]

Ok, the rationale of allowing "any complete slicing" (as opposed to partial slicing, which requires more allocations than just the unavoidable unit vector and the result) sounds convincing to me.

As I said, it's borderline and a matter of taste: from the point of view of matrix construction it's clear that this is expensive (and could deserve a warning), from the point of view of slicing it is "expectedly expensive" (and could therefore be allowed). I'll see if I need to adjust the announcement. Perhaps, there should also be a section in the types.md docs, that repeats considerations and lays out the rationale of what is allowed and what not.

[dkarrasch]

Uff, I think I finally understood the macro stuff! 🎉

What shall we do then? Allow slicing in general without interference of macros and allow scalar indexing/partial slicing upon a macro signal, or allow any slicing/indexing only upon macro signal? Having two macros, one for slicing and one for indexing and partial slicing is perhaps overkill.

[oschulz]

What shall we do then? Allow slicing in general without interference of macros

From a user perspective, I'd be happy with A[:,j] and A[:,:] allowed, A[i,:] allowed if adjoint exits, and everything else resulting in an error with a helpful "don't do this, it's slow" error message. Anything beyond (macro magic to allow stuff explicitly, etc.) is fine too, of course.

[dkarrasch]

Ok, then you should be happy with the current state, I hope. We could release this as is and then take some time to work out how to allow things that are currently forbidden. That would be "a new feature" and non-breaking, as opposed to first allowing something and later require a macro to do that. I'll wait for @JeffFessler's thumb.

[oschulz]

Sounds great - again, thanks a lot!

[JeffFessler]

I noticed that codecov complains about these two lines. I wonder if they should be commented out because we are not supporting A[1] and A[end] at this point in time, right?

[JeffFessler]

should these be commented out since A[I,J] is not supported?
the codecov warning is what caught my eye here.

[dkarrasch]

The tests were not variable enough. This method is needed in general.

[JeffFessler]

I'd be happy with A[:,j] and A[:,:] allowed, A[i,:] allowed if adjoint exits

this is also fine with me. thumbs up other than a couple questions i asked about commenting out possibly unused code which i leave to you to decide.

[oschulz]

Thanks again!