Make LinearMap a subtype of AbstractMatrix to increase code compatibility?

[oschulz]

My apologies in case this has been discussed before. @dkarrasch would you consider making LinearMap{T} a subtype of AbstractMatrix{T}?

I've really grow to love the flexibility of LinearMaps and the ability to use them as drop-in replacements for matrices in many applications. But sometimes that's not easy because generic code in packages frequently requires an AbstractMatrix as the argument type - a LinearMap would often work well since a lot of code only uses high-level linear algebra operations, but LinearMap and AbstractMatrix just have no common supertype that typical package code would dispatch on.

To me, a LinearMap feels like a matrix in almost every way already - is has a size, element type, etc. What it currently doesn't have is getindex - but I think it could. A[:,:] and A[:,j] would be efficient, and if A is transposable then A[i,:] would be efficient as well. Even views would be somewhat efficient I guess, just A[i,j] would be horribly inefficient of course. But the same is true for Julia GPU arrays like CuArray - one could adopt the same warning/error mechanism implemented there that alerts users that serial/element-wise access is inefficient (or even prevents it).

[oschulz]

It would be neat to support (e.g.) A \ b as well, but that would require depending on IterativeSolvers.jl (would increase the load time of LinearMaps from about 110 ms to about 140 ms) or Krylov.jl (seems a bit more lightweight).

Alternatively, we could have a wrapper AsMatrix(::LinearMap) in a separate package (call it ArrayMaps.jl or so?). That wrapper could then also take solver type/parameters as an optional argument (a good default setting could be chosen automatically since LinearMap already provides information on whether it's pos-def and so on).

(Example use case: we do solve A * x == b in MGVI.jl a lot (via conjugate gradient), we just can't write it nice as A \ b yet :-) )

[oschulz]

It would be neat to support (e.g.) A \ b as well [...] LinearMap already provides information on whether it's pos-def and so on

I asked @ChrisRackauckas for advice, good defaults might me Krylov.gc if posdef, Krylov.symlq (or Krylov.minres?) if symmetic, and Krylov.gmres otherwise. Krylov.jl is very lightweight and dependency-free - could it maybe be an acceptable dependency for LinearMaps to take on?

julia> import DensityInterface # To get some Pkg init stuff out of the way

julia> @time_imports using LinearMaps, Krylov
    104.9 ms  LinearMaps
     10.3 ms  Krylov

[JeffFessler]

This was discussed in #38. Nice to see another vote in favor of AbstractMatrix though.
See https://github.com/JeffFessler/LinearMapsAA.jl

[oschulz]

Thanks for the link @JeffFessler , I hadn't been aware of LinearMapsAA and I had overlooked #38.

I do get the arguments in #38 regarding getindex - but an ineffecient getindex hasn't stopped GPU arrays from being very successful. There is a lot of code that asks for an AbstractArray or AbstractMatrix but would work efficiently with LinearMaps since it just uses high-level linear algebra. I feel that would justify making LinearMap a subtype of AbstractMatrix even if arbirary getindex wouldn't be efficient (or single-element access not even implemented). In the end, that would probably also be more lightweight and efficient than having to add an additional wrapper layer.

[oschulz]

There's also the issue that code that want's to support both AbstractMatrix and LinearMap currently has to depend on LinearMaps.jl, which has a load time of about 100 ms. Not super-heavy, but also not lightweight - many lighter packages will want to avoid it.

[Jutho]

@oschulz , is there any reason to "depend on" LinearMaps.jl for compatibility? Why not just be generic and leave argument types unspecified; any user input that can multiply with a vector and produce a new vector would be accepted? That was at least the advice that was given to me by Steven Johnson when I started using Julia in 2014 and started working on this package.

The choice for not making LinearMap a subtype of AbstractMatrix back then was because an efficient (i.e. O(1) ) getindex seemed to be part of the AbstractArray interface. CuArray does indeed seem to be a somewhat borderline case. As to making LinearMap a subtype of AbstractMatrix, I leave it to others, in particular @dkarrasch , to decide whether that decision needs to be revised.

It seems the lack of formal interface specifications, and the potential issues that arise from that, start to pop up in various places in the Julia ecosystem.

[oschulz]

@oschulz , is there any reason to "depend on" LinearMaps.jl for compatibility? Why not just be generic and leave argument types unspecified

But that's basically saying "don't use multiple dispatch". Code may need to distinguish between "array-like" object and other types, e.g. scalars, to implement optimized methods. The AbstractArray supertype (plus AbstractMatrix, etc.) has been (and I think made) Julia very successful. Not being able to express that I need something that behaves in a matrix/operator-like fashion is limiting. And one couldn't ensure match of number of dimension, or dispatch on element type, or ...

It seems the lack of formal interface specifications, and the potential issues that arise from that, start to pop up in various places in the Julia ecosystem.

I couldn't agree more - and it leads to a web of complicated (semver-unsafe) @requires or lightweight packages being made heavy by large dependencies, etc. In some cases, using pure interfaces without supertypes can be a solution, but I don't think that's what we should to in general, because we'd lose all the power of multiple dispatch.

Since Julia doesn't have multiple dispatch (directly), we can't really have an "operator-like" supertype for AbstractMatrix{<:Number}, and using Holy-traits all over linear-algebra code would be messy too. That why I'd argue that making LinearMap and AbstractMatrix would be a clean solution. I would say that (again looking at GPU array) there's a community consensus that single-element getindex of an array doesn't have to be efficient (with appropriate guards in place).

[Jutho]

But that's basically saying "don't use multiple dispatch". Code may need to distinguish between "array-like" object and other types, e.g. scalars, to implement optimized methods. The AbstractArray supertype (plus AbstractMatrix, etc.) has been (and I think made) Julia very successful. Not being able to express that I need something that behaves in a matrix/operator-like fashion is limiting. And one couldn't ensure match of number of dimension, or dispatch on element type, or ...

In an ideal world (which I realise is not what we always live in), that optimisation should be done at the level of the methods that implement "linear map (of any kind)" * vector, no? By scalar I assume you mean, scalar times the identity matrix? I do assume you don't want to have completely different behaviour if the argument is a scalar instead of a matrix? Cause that would be abusing multiple dispatch; completely different behaviour would imply different methods.

Note that LinearMap objects do listen to size and eltype. It really seems that, what is missing, is a proper interface for what it means to be a linear map, like there is the iteration interface (which is also not associated with any particular type hierarchy).

It is not that I am radically against making LinearMap <: AbstractMatrix, it is just that it seems that whatever methods you have in mind, may also work with other types that behave as linear maps. So it seems to me (but of course I don't know your specific use case) that it should be possible to have a number of specialised methods, where the argument is a Number or a Matrix, and then a catch all method with unspecified argument type, that imposes a minimum of assumptions on this argument, namely something like size, eltype, * with vector.

[oschulz]

By scalar I assume you mean, scalar times the identity matrix?

For example, but also more general let's say you have something like

mcmc_sample(target, metric::Real} =...
mcmc_sample(target, metric::AbstractMatrix{<:Real}} = ...
mcmc_sample(target, metric::MLAutoAdaptMetic} = ... # very different implementation

A LinearMap would work as a metric, but I can't pass it. So I could define

mcmc_sample(target, metric} = ...

for the case that metric is something like an array or linear map. But it'll be hard to see in the code. And also, what if there's the same situation in a hypothetical "auto-adapt-metic" ecosystem - several different ones with the same duck-typed interface but no common supertype? For that I'd need mcmc_sample(target, metric) without type-restriction as well. But the implementation is very different from the operator/matrix-like case (let's just assume it is). So then I start with if-method-defined, etc. - essentially ending up with Python and loss of multiple-dispatch and composability.

It is not that I am radically against making LinearMap <: AbstractMatrix, it is just that it seems that whatever methods you have in mind, may also work with other types that behave as linear maps.

Well, just on the purely practical side I guess many algorithms in packages (and possibly std libs as well) that could efficiently utilize LinearMaps, but we won't be able to make them drop the ::AbstractMatrix in their methods signatures. :-)

Also, just speaking for myself, I do like restricting method arguments - types carry semantic meaning, and can make code more readable and wrong types get caught earlier. And multiple-dispatch needs good type hierarchies.

[oschulz]

it is just that it seems that whatever methods you have in mind, may also work with other types that behave as linear maps.

True, there's several "linear operator/map" kinds around ... we have LinearMaps.jl, LinearOperators.jl (very similar but different type - maybe some synergy potential there), FFTW.jl plans support at least part of the "linear operator interface" (CC @stevengj) and maybe there's more as well.

It would be nice to have an AbstractLinearOperator/AbstractLinearMap super-type in a lightweight package, as well as an interface definition. Of course duck typing can be very successful, but you can't dispatch on it (easily, without Holy traits). Of course in some situations, one needs traits/duck-typing because the natural/primary super-type of types lies in another domain than the one an interface want's to address - Julia doesn't have multiple inheritance, after all. But when the "primary domain" of a type is clear I feel that a common type hierarchy is pretty much always beneficial, since then one doesn't need wrapper/adaptor types - the benefit of multiple dispatch.

[oschulz]

CC @dpo, @abelsiqueira and @geoffroyleconte (LinearOperators.jl).

[JeffFessler]

When I first brought up #38, I also felt pretty strongly that a LinearMap should be a subtype of AbstractMatrix, which led to LinearMapsAA. Over time, however, I found that what I really need in most of my work is more general linear operators that map from arrays to arrays, rather than from vectors to vectors. For example, the Radon transform in tomography maps 2D images into 2D sinograms. So in LinearMapsAA there are constructors for another type (LinearMapAO where the O is for operator) and that type cannot be a subtype of AbstractMatrix because it is design to be multiplied with arrays, not with vectors. Another example is the 2D FFT that maps a 2d array into another 2d array. The package https://github.com/hakkelt/FunctionOperators.jl also supports mapping arrays to arrays.

I agree that it could be beneficial to have a AbstractLinearMap super-type.