Support for lazily sized maps?

[marius311]

I've been toying with replacing some of my hand-written LinearMap-like code with this package, which is much nicer than what I have. But one roadblock is that to be really useful, I'd need support for "lazily sized" or "unsized" function maps. Basically, a map which can implicilty take on any size it needs to depending on the vector its being applied to. The FFT is one example. While I can write,

A = LinearMap(fft, 4, 4, ismutating=false)
A * rand(4)

if I do

A * rand(5)

its a dimension error, even though the code would work if it just didn't do the dimension check

I'm curious whether there's any interest in adding something to let this work? Quite possibly it breaks some other more complex things LinearMaps does? One interface which seems reasonable to me would be to let the user pass : for the dimension (Any might be another option), like

LinearMap(fft, :, :, ismutating=false)

and then the dimension check code could be:

A.N === (:) || length(x) == A.N || throw(DimensionMismatch())

This would clearly only be possible for ismutating=false maps.

[Jutho]

The original motivation for this package was that eigs from Arpack.jl (back in the days when this was still included in Julia Base) required a "matrix-like" object for its linear operator of which it wants to compute eigenvalues. In particular, that object needs to listen to size and object * vector. That's why linear map objects from this package have a fixed size. I am not sure but I think all the stuff in IterativeSolvers.jl is also like this.

If really all you have is a function that you would like to use in a Krylov algorithm, you could try KrylovKit.jl, which even doesn't care about the size of the vectors that you combine it with.

[marius311]

My use-case isn't so much that I already have a function, but that its really convenient to use this package's composition to build one, before feeding it into various algorithms, including Kyrlov ones. As an example, say I want to Wiener filter some data where my signal is diagonal in Fourier space. I'm envisioning doing something like:

signal_spectrum = ... # some vector
noise_variance = ... # some vector

S = LinearMap(Diagonal(signal_spectrum))
N = LinearMap(Diagonal(noise_variance))
F = LinearMap(fft, ifft, :, :, ismutating=false)

(F' * S * F) * conjugate_gradient(F' * S * F + N, data)  # compute Wiener filter

and it would be really nice not to have to specify the size of F so I can reuse it for all problem sizes (this is already what my hand-rolled thing does, just with code that's way less nice than LinearMaps). Maybe this example makes it seem too easy to just suck it up and specify the size for F each time, but in my real case I have a lot more of these lazily sized operators floating around to where the tradeoff of having to specify their sizes wouldn't really be worth it for me to switch to LinearMaps.

Anyway, just want to give a flavor of what I'm doing, I can understand if ultimately you don't think its a good fit for this package.

[Jutho]

So do the vectors signal_spectrum and noise_variance have a fixed length? Why can't you just use F = LinearMap(fft, ifft, length(signal_spectrum), ismutating = false). If not, I guess I see two options:

1. all the linear maps in your example have a flexible size. I don't see how you would be able to use the resulting linear maps in any of the typical packages that you want to combine it with, except for KrylovKit

2. or in your example only F has a flexible size, and it acquires a fixed size when combining it with other linear maps which specify the size. That would be possible, theoretically, but hard to implement in the current implementation, and it would still disallow you from using F alone in e.g. IterativeSolvers.jl or Arpack.jl.

Furthermore, in the case 2, I don't see why not just defining F with fixed size using length(signal_spectrum) would not be possible?

[marius311]

Yea, its exactly case 2, and the answer to your quesion is what I was referring to with

Maybe this example makes it seem too easy to just suck it up and specify the size for F each time, but in my real case I have a lot more of these lazily sized operators floating around to where the tradeoff of having to specify their sizes wouldn't really be worth it for me to switch to LinearMaps.

Some other examples of operators I have that are lazily sized are a spatial gradient operator (bound to ∇), a resolution upgrade/downgrade operator, or a low/high pass operator that filters out modes above/below a certain Fourier wave number. The code to construct any of these does not need to know the size of the input vector at operator creation time. Again, its not impossible to rewrite everything to pass around explicit sizes, but probably not worth it in my case.

Maybe I fork this and play around with it for a while, and come back if I have a more definite proposal. I'm not super familiar with this code but on quick glance it kind of seems like loosening FunctionMap's requirement that M/N::Int to M/N::Union{Int,Colon} (in a type stable way) and modifying the bounds checks would get you very far. Another option could be if there's a halfway point where this package allows : as a dimension specifier but doesn't do anything with it (which could be a non-breaking change), but then this allows an optional separate package which makes the lazily-sized machinery work.

[JeffFessler]

My own work is typically in 2D and 3D so fft has even more generality there.
That is why https://github.com/JeffFessler/LinearMapsAA.jl supports a version that maps arrays to arrays, not just vectors to vectors. I am intrigued by this lazy sizing idea, but generalizing it to arrays would take some thought because if you do
Afft * rand(8,7) do we want the 2D FFT of the RHS, or the 1D FFT of each column of the RHS? Specifying the input dims avoids this ambiguity. Perhaps LinearMapAA(fft, (:,:), (:,:)) could convey that this is a lazy 2D FFT. @marius311 you mention terms that sound like you might also need 2D operations?

[JeffFessler]

Also, operations like up-sampling are non-square so I think there would need to be a way to describe the "aspect ratio" for this to all work with other operations like cat.

[marius311]

My work is indeed 2D (maps of the Cosmic Microwave Background). I had seen LinearMapsAA but missed the arrays to arrays part. In what I do, I handle this with a custom array type which is pretty similar to ComponentArrays, in that it splays out a Matrix into a Vector representation.

The example from you link is one I should have mentioned here too, lazily sized maps is exactly what I already is in Base.