JuliaLang / julia

The Julia Programming Language
https://julialang.org/
MIT License
45.93k stars 5.49k forks source link

permutedims(::SparseMatrixCSC) calls a fallback method #28324

Closed jlapeyre closed 5 years ago

jlapeyre commented 6 years ago

The documentation recommends permutedims to get a material transpose. There is an efficient function, but it is not called:

julia> size(S)
(20, 10)

julia> typeof(S)
SparseMatrixCSC{Float64,Int64}

julia> @btime permutedims($S);
  2.860 μs (12 allocations: 2.38 KiB)

julia> @btime copy(Transpose($S));
  188.709 ns (5 allocations: 720 bytes)

julia> copy(Transpose(S)) == permutedims(S)
true
andreasnoack commented 6 years ago

We'll need a permutedims methods for SparseMatrixCSC that just calls copy(transpose(A)). Would be great if you could prepare a PR.

jlapeyre commented 6 years ago

Would be great if you could prepare a PR.

OK. In case someone cares to comment first, these are the relevant lines in sparsematrix.jl

adjoint(A::SparseMatrixCSC) = Adjoint(A)
transpose(A::SparseMatrixCSC) = Transpose(A)
Base.copy(A::Adjoint{<:Any,<:SparseMatrixCSC}) = ftranspose(A.parent, conj)
Base.copy(A::Transpose{<:Any,<:SparseMatrixCSC}) = ftranspose(A.parent, identity)

The following should provide a minimal fix

Base.permutedims(A::SparseMatrixCSC) = ftranspose(A.parent, identity)
function Base.permutedims(A::SparseMatrixCSC, (a,b))
    (a, b) == (2, 1) && return ftranspose(A.parent, identity)
    (a, b) == (1, 2) && return copy(A)
    throw(ArgumentError("no valid permutation of dimensions"))
end
andreasnoack commented 6 years ago

You are right that permutedims is only supposed to act on the top-level so my proposed solution wasn't correct.

Using copy(transpose(M)) to get a material transpose is so obviously wrong...This use is semantically orthogonal to the other uses of copy.

This has indeed been discussed. Maybe @Sacha0 remembers which are the relevant issues. Would you prefer to get a Transpose when using copy? Would you copy the underlying array or should it just be a noop? A similar example where it might be more obvious that you don't want copy to return the same type is for SubArray where it would be unfortunate if copy(view(randn(1000,1000), 1:2, 1:2)) copied the underlying arrays and it is redundant to wrap the copied 2x2 matrix in a SubArray.

jlapeyre commented 6 years ago

it would be unfortunate if copy(view(randn(1000,1000), 1:2, 1:2)) copied the underlying arrays and it is redundant to wrap the copied 2x2 matrix in a SubArray

If it comes down to a choice between preserving the current consensus for copy(transpose(.... and disasters like your example (or the absurdity of copying just a wrapper), then clearly the former is better. It may be non-orthogonal design (it certainly was not "discoverable" for me), but it is performant.

Would you prefer to get a Transpose when using copy? Would you copy the underlying array or should it just be a noop?

That's a good question. But, maybe better questions are what do you want and what is the best interface for providing it ? Currently, we have

julia> v = view(randn(1000, 1000), 1:2, 1:2);

julia> copy(v)
2×2 Array{Float64,2}:
 -0.981013  1.06454
 -1.69617   1.73077

julia> collect(v)
2×2 Array{Float64,2}:
 -0.981013  1.06454
 -1.69617   1.73077

The same holds for collect(transpose(m)). It seems to me (again, I'd like to see more counter examples) that recommending collect for materializing the view is better precisely because it is less ambiguous. So, the current behavior of copy might be correct for avoiding unintended consequences. But, that doesn't mean it has to be the recommended interface for materializing lazy views. If the behavior of collect above holds more or less generally, then it shouldn't be a huge effort to switch to recommending it. Has using collect this way been discussed ?

andreasnoack commented 6 years ago

Has using collect this way been discussed ?

Yes and for some reason it has fallen in dismay. I'm pretty neutral wrt collect but the general view seems to be that it should go away and in favor of using constructors. Most often that would be Array. However, in this case, and possibly other cases, something generic is needed since the result should depend on whatever is being wrapped by Transpose. Again, @Sacha0 has a good memory and will probably be able to fill in the missing details here.

jlapeyre commented 6 years ago

Ok thanks for the tip. Here are some issues to start with: JuliaLang/julia#16029, JuliaLang/julia#24595, JuliaLang/LinearAlgebra.jl#231. Looks like it started with more or less "collect can always be replaced with Vector". But it turns out it's more complicated than that.

Sacha0 commented 6 years ago

Regarding the semantics of copy, collect, materialize and friends, I

punt

to the inimitable @mbauman. Best!

jlapeyre commented 6 years ago

I punt

I'm preparing a what I hope is a non-controversial PR :)

Sacha0 commented 6 years ago

I'm preparing a what I hope is a non-controversial PR :)

Famous last words 😉.

mbauman commented 6 years ago

Yes, copy vs. collect vs. materialize goes down a big rabbit hole very quickly. If you're curious about reading more, JuliaLang/julia#16029 is a good starting point for the latest thoughts. From there leads a labyrinth of linked issues.

For now, we're living with the fact that copy plays two roles: it copies elements into a new outer structure and in doing so it often "materializes" laziness. This is something I know both Sacha and I have wanted to address, but it's a big task that's been on hold as we've finalized 0.7/1.0.

ViralBShah commented 5 years ago

This appears resolved to me. I get the same performance for both operations reported above.