recognize Adjoint/Transpose of sparse array as sparse

[dkarrasch]

Closes #34253.

[ViralBShah]

It would be good to add a test for the sparse matrix adjoint case too, which has the same bug, which I believe this PR should fix.

[dkarrasch]

It's ten lines above these. ;-)

[oxinabox]

Consider also:

julia> adjoint(transpose((1 + im) .*sprandn(5,5, 0.5)))
5×5 Adjoint{Complex{Float64},Transpose{Complex{Float64},SparseMatrixCSC{Complex{Float64},Int64}}}:
  0.103342-0.103342im        0.0-0.0im       -0.206055+0.206055im  -0.559357+0.559357im     2.0501-2.0501im
 -0.558443+0.558443im   -1.16483+1.16483im   -0.407952+0.407952im        0.0-0.0im             0.0-0.0im
  0.564176-0.564176im        0.0-0.0im             0.0-0.0im             0.0-0.0im        0.465629-0.465629im
  0.117705-0.117705im  -0.545153+0.545153im        0.0-0.0im             0.0-0.0im       -0.469461+0.469461im
       0.0-0.0im       -0.410911+0.410911im  -0.364878+0.364878im        0.0-0.0im             0.0-0.0im

You may want to do:

issparse(S::LinearAlgebra.Adjoint) where T = issparse(parent(S))
issparse(S::LinearAlgebra.Transpose) = issparse(parent(S))

However, this may trigger a runtime call to parent?
It shouldn't since it only needs the type that is returned by parent
which should be trivially inferred.

If it isn't might be good to make these proper traits that take only the type,
not an instance.
and then just have a outer issparse(x) = issparse(typeof(x))

---

Though perhaps
Adjoint{T ,Transpose{T,SparseMatrixCSC{T,Int64}}}
actually should not be considered sparse, for purposes of issparse

[dkarrasch]

Though perhaps
Adjoint{T ,Transpose{T,SparseMatrixCSC{T,Int64}}}
actually should not be considered sparse, for purposes of issparse

What do you mean? I liked the trait idea very much, actually. That would also resolve stuff like this.

julia> A = UpperTriangular(sprand(10,10,0.2))';

julia> issparse(A)
false

Should we merge this as is (because it closes another concrete issue) and open a new issue to discuss how many layers of wrapping should be considered as sparsity-preserving?

[ViralBShah]

I am in favour of merging this one and opening a new one for the larger sparsity preserving discussion.

[oxinabox]

Though perhaps Adjoint{T ,Transpose{T,SparseMatrixCSC{T,Int64}}} actually should not be considered sparse, for purposes of issparse

What do you mean?

Well it depends what one is promising to do by returning issparse(x)=true.
Like do you promise that a certain set of functions exists, and has a certain performance characteristic?
Seems not, 🤷‍♂ might as well say it's sparse, cos conceptually it is

[oxinabox]

I just checked: the parent one constant folds.
I feel like it is a strict improvement over the current PR,
and is an easy change to make.

harder versions like doing this for all arrays that have parent defined, another PR.
Or changing it all to traits, another PR.

[dkarrasch]

Would

issparse(S::LinearAlgebra.Symmetric) = issparse(S.data)

and alike be also acceptable? That would make issparse almost a trait, right? Is there any additional runtime cost to be expected?

[oxinabox]

Not that issparse(S::LinearAlgebra.Symmetric) = issparse(S.data)
is the same as issparse(S::LinearAlgebra.Symmetric) = issparse(parent(S))

AFAICT, there is a convention (which may not have ever been written down).
that wrapper arrays store the inner array in the .data field.
And thus that parent(x) = x.data.

Is there any additional runtime cost to be expected?

Not in this case.

julia> SparseArrays.issparse(S::LinearAlgebra.Symmetric) = issparse(parent(S))

julia> x = Symmetric(sprand(5, 5, 0.5));

julia> @code_typed issparse(x)
CodeInfo(
1 ─     return true
) => Bool

[dkarrasch]

Oh, I see, so I could extend the parent approach to all these:

issparse(S::LinearAlgebra.Symmetric{<:Any,<:AbstractSparseMatrix}) = true
issparse(S::LinearAlgebra.Hermitian{<:Any,<:AbstractSparseMatrix}) = true
issparse(S::LinearAlgebra.LowerTriangular{<:Any,<:AbstractSparseMatrix}) = true
issparse(S::LinearAlgebra.UnitLowerTriangular{<:Any,<:AbstractSparseMatrix}) = true
issparse(S::LinearAlgebra.UpperTriangular{<:Any,<:AbstractSparseMatrix}) = true
issparse(S::LinearAlgebra.UnitUpperTriangular{<:Any,<:AbstractSparseMatrix}) = true

?

EDIT: The parent approach seems to successfully constant fold through two layers, and likely more.

[oxinabox]

Yep, and most of:

julia> methods(parent)
# 10 methods for generic function "parent":
[1] parent(V::SubArray) in Base at subarray.jl:79
[2] parent(A::Base.ReshapedArray) in Base at reshapedarray.jl:209
[3] parent(a::Base.ReinterpretArray) in Base at reinterpretarray.jl:97
[4] parent(A::PermutedDimsArray) in Base.PermutedDimsArrays at permuteddimsarray.jl:48
[5] parent(A::Union{Adjoint{T,S}, Transpose{T,S}} where S where T) in LinearAlgebra at /usr/local/src/julia/julia-master/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/adjtrans.jl:200
[6] parent(A::LinearAlgebra.AbstractTriangular) in LinearAlgebra at /usr/local/src/julia/julia-master/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/triangular.jl:163
[7] parent(A::Union{Hermitian{T,S}, Symmetric{T,S}} where S where T) in LinearAlgebra at /usr/local/src/julia/julia-master/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/symmetric.jl:270
[8] parent(D::Diagonal) in LinearAlgebra at /usr/local/src/julia/julia-master/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/diagonal.jl:112
[9] parent(H::UpperHessenberg) in LinearAlgebra at /usr/local/src/julia/julia-master/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/hessenberg.jl:57
[10] parent(a::AbstractArray) in Base at abstractarray.jl:1136

Not the last though, that will infinite loop.
Since non-wrapper arrays default to being own parent.

But on that basis idea, get all the wrapper arrays:

@inline function SparseArrays.issparse(x::AbstractArray)
    p = parent(x)
    return p===x ? false : issparse(p)
end

Which also still constant folds away.
Need to check @code_llvm for that though, as it needs LLVM's optimization passes not just julia's

julia> z = adjoint(transpose((1 + im) .*sprandn(5,5, 0.5)));

julia> @code_llvm issparse(z)

;  @ REPL[18]:2 within `issparse'
define i8 @julia_issparse_17287(%jl_value_t addrspace(10)* nonnull align 8 dereferenceable(8)) {
top:
;  @ REPL[18]:4 within `issparse'
  ret i8 1
}

julia> q = adjoint(transpose((1 + im) .*randn(5,5)));

julia> @code_llvm issparse(q)

;  @ REPL[18]:2 within `issparse'
define i8 @julia_issparse_17303(%jl_value_t addrspace(10)* nonnull align 8 dereferenceable(8)) {
top:
;  @ REPL[18]:4 within `issparse'
  ret i8 0
}

[dkarrasch]

@oxinabox So what's the conclusion? I'm not familiar with LLVM optimization and don't know how to read @code_llvm, so do you think I should include your proposal?

[oxinabox]

@ViralBShah was right, this is too much.
You make this PR as is and I will make a follow up.