How to handle Linear operators / functionals?

[dlfivefifty]

At the moment we have the following setup for linear operators / functionals that are not built out of elementary operations, where f isa AbstractQuasiVector and T isa AbstractQuasiMatrix:

1. Derivative(x) * f and Derivative(x) * T
2. sum(f) and sum(T; dims=1)
   As we consider other linear operators such as cumsum it might be a good idea to see how to make this consistent

Version (2) is in some sense "canonical" as it exists in Base. So the question is what to do with Version (1). We could have something like

struct Linear{T,Axes,F} <: AbstractQuasiMatrix{T}
    f::F
    axes::Axes
end

*(L::Linear, v::AbstractQuasiVector) = L.f(v)
*(L::Linear, v::AbstractQuasiMatrix) = L.f(v; dims=1)
axes(L::Linear) = L.axes

const Diff{T,Axis} = Linear{T, Axis,typeof(diff)} # Replaces Derivative
const Cumsum{T,Axis} = Linear{T, Axis,typeof(cumsum)}
const Sum{T,Axis} = Linear{T,Axis,typeof(sum)}

Diff(x) = Linear{Float64}(diff, (x,x))
Cumsum(x) = Linear{Float64}(cumsum, (x,x))
Sum(x) = Linear{Float64}(diff, (Base.OneTo(1),x))

How this relates to axes-free alternatives is another question: #22

[dlfivefifty]

The reason we like having Derivative is for making operators like

(ε*D^2 + x .* D + I) * T

In theory an option is to no longer support this and force users to do

ε*diff(diff(T;dims=1); dims=1)  + x .* diff(T; dims=1) + T

[jagot]

I would be very unhappy with version (2), not least because I would have to rewrite a lot of code, but also because I feel it is much uglier. I guess supporting this as an alternative is fine.

How would you apply the transpose of an operator this way?

[dlfivefifty]

Right that wasn't really a serious proposition. But we should support at minimum both

[jagot]

My use case is usually R'*A*R where R isa Basis and A is some kind of linear operator, so I'm mostly interested in generating matrix representations of linear operators. I have not yet worked with storing linear operators lazily and applying them on-the-fly; I'm assuming that in the example above T = R*c for some vector c of expansion coefficients.