size() returns Tuple of Tuples instead of Tuple

[Benoit9]

I am using the Conv operator as the first argument matrix-like operator in the Krylov.jl optimization functions.

It fails because the size() function applied to a Conv object does not return the expected tuple. Instead of

(m, n)

The Krylov.jl optimizer gets:

((m,), (n,))

If I override Base.size, it all works beautifully. Shouldn't the size operator conform to what is expected of a Matrix-like operator? Am I doing something wrong?

Edit: in order to silence a warning, I also had to override eltype() to return Float64 (in my case) instead of Any.

[nantonel]

Hi, yes as far as I remember in this package size always returns tuples of tuples. This is because some operators not only work with matrices but with tensors too, e.g. size can be ((m1,m2,...),(n1,n2,...)).

Instead of overriding the functions I suspect it would be easier to create a function in Krylov.jlwith dispatch <: AbstractOperators where you extract (n,m) from ((m,), (n,)) and use codomainType instead of eltype.

Alternatively you could modify the function size(A::AbstractOperator) to splat ((m,), (n,)) to (n,m) if:

all( [length(s) for s in size(A)] .== 1)

but I think that might break some code in here, so it might be a time consuming change.

[Benoit9]

Thank you for the quick and detailed answer!

Is it usual practice to have size() return a Tuple of Tuples in the case of tensors? I am not familiar enough with the Julia ecosystem to tell. Should I open an issue on Krylov.jl instead?

I was thrilled to be able to just plug in AbstractOperators.Conv to a Krylov.jl solver. The override of size() is just a minor wrinkle.

[nantonel]

Happy you're enjoying this package!

Is it usual practice to have size() return a Tuple of Tuples in the case of tensors?

The thing is that in Julia a matrix size is:

julia> size(randn(10,10))
(10, 10)

and a tensor is:

julia> size(randn(2,3,4))
(2, 3, 4)

In this package we defined the size of the abstract operator as:

Returns the size of an AbstractOperator. Type size(A,1) for the size of the codomain and size(A,2) for the size of the codomain.

I have no clue if the devs of Krylov.jl are interested in supporting this package but you can certainly give it a shot!

[hakkelt]

It's also a solution to wrap the AbstractOperator with a LinearOperator. For example:

op = Conv(rand(10), rand(3))
dft = LinearOperator(eltype(op), prod(size(op, 1)), prod(size(op, 2)), false, false,
                     (x, y) -> mul!(x, op, y), 
                     nothing,  # will be inferred
                     (x, y) -> mul!(x, op', y))

I might add a function that does this wrapping automatically sometime soon as a package extension if you don't have an objection, @nantonel.