Julia: Inconsistency in descending `transpose` between sparse and dense arrays

[jlapeyre]

using LinearAlgebra
using SparseArrays
A = [1 0 ; 1 0]

julia> M = reshape([A, A, A, A], (2,2))
2×2 Array{Array{Int64,2},2}:
 [1 0; 1 0]  [1 0; 1 0]
 [1 0; 1 0]  [1 0; 1 0]

julia> copy(transpose(M))
2×2 Array{Array{Int64,2},2}:
 [1 1; 0 0]  [1 1; 0 0]
 [1 1; 0 0]  [1 1; 0 0]

julia> Array(copy(transpose(sparse(M))))
2×2 Array{Array{Int64,2},2}:
 [1 0; 1 0]  [1 0; 1 0]
 [1 0; 1 0]  [1 0; 1 0]

Ideally, I expect sparse and dense objects to represent the same thing, but to employ different representations internally. I convert between them only to take advantage of differences in algorithmic efficiency.

I guess there is no perfect solution to this kind of thing. Computer-language types and mathematical-object types (sets) are far from congruent. E.g., in math, it is common to interpret objects in one section of a manuscript or book in a particular way without worrying about breaking examples in another section. Even the most flexible computer-language type system is more rigid.

[andreasnoack]

It's a bug. Maybe easier to see in

julia> transpose(sparse(M))[1,1]
2×2 Transpose{Int64,Array{Int64,2}}:
 1  1
 0  0

julia> copy(transpose(sparse(M)))[1,1]
2×2 Array{Int64,2}:
 1  0
 1  0

I guess copy(::Transpose/Adjoint) needs to be adjusted somehow.

[andreasnoack]

Indeed the two definitions https://github.com/JuliaLang/julia/blob/81ee3ae0cad505b7e82115820a6e5a7116019137/stdlib/SparseArrays/src/sparsematrix.jl#L870-L871 are only correct for scalar element types.

[fredrikekre]

IIRC the assumption of scalar elements are used quite a lot in the sparse code.

[jlapeyre]

Not sure how robust this is, but

julia> Base.copy(A::Transpose{<:Any,<:SparseMatrixCSC}) = 
     SparseArrays.ftranspose(A.parent,  x -> copy(transpose(x)));

julia> S = sparse(M);

julia> copy(transpose(S))[1,1]
2×2 Array{Int64,2}:
 1  1
 0  0