Should hcat of dense and sparse matrix be sparse?

[pkofod]

I wanted to hcat a dense vector B and a sparse matrix A. I expected

n = 5
A = speye(n)
B = ones(n)
C = [B A]

To return a sparse matrix C, but instead A is converted to it's dense representation, and this means C is obviously also dense. I then did

C = sparse([B A])

But as @tavert pointed out over at reddit, it is better to.

C = [sparse(B) A]

My question is, it is intended that users should do the latter, or is the current behaviour a bug or at least unintended?

[mlubin]

Definitely not a bug, in some cases it's the reasonable thing to do. It really depends on the:

1. Density of the "dense" block
2. Relative dimensions
3. What you're going to do with the matrix afterwards.

None of the three above are properties of the types themselves, so to do something intelligent (if that's possible at all) you would lose type stability. I vote for the simple and predictable behavior: sparse + dense = dense.

[tkelman]

Matlab and Octave return sparse in this case. Would have to check scipy.

edit: Julia's currently going to a generic fallback here:

julia> @which hcat(ones(5), speye(5,5))
hcat{T}(A::Union{AbstractArray{T,1},AbstractArray{T,2}}...) at abstractarray.jl:742

[ViralBShah]

While there are different things for different situations, I think the default should be to preserve sparsity in cases of concatenation.

[tbreloff]

I tend to agree with Viral. If you're using sparse matrices then you
already have your eye toward storage, so that probably makes the most
sense. Returning dense or sparse depending on matrix features is just
asking for trouble.

On Mon, Sep 14, 2015 at 11:00 PM, Viral B. Shah [email protected]
wrote:

While there are different things for different situations, I think the
default should be to preserve sparsity in cases of concatenation.

—
Reply to this email directly or view it on GitHub
#13130 (comment).

[ViralBShah]

This requires adding a new definition that calls sparse on all AbstractArray inputs, so long as one of the inputs is sparse.

[pkofod]

Out of curiosity (I would love to do it myself, but I'm not sure how it would work), how would you do that with varargs? I see something like function hcat(X::SparseMatrixCSC...) at https://github.com/JuliaLang/julia/blob/master/base/sparse/sparsematrix.jl#L2478 , but how do you identify that at least one object in a vararg is of type SparseMatrixCSC? Check all hcats of AbstractArrays?

[ViralBShah]

Union{Matrix,SparseMatrixCSC} ?

[pkofod]

Well, of course :) I'll have a go at it.

[pkofod]

Well, I guess I do not totally understand your suggestion @ViralBShah . If I add a definition of hcat and/or vcat as function vcat(X::Union{Matrix, SparseMatrixCSC}...) it will make all concatenation of matrices (dense or sparse) sparse. Maybe this is not what you meant? Union gives "at least one of Matrix or at least one of SparseMatrixCSC" but I need "at least one of Matrix and at least one of SparseMatrixCSC".

[tkelman]

I don't know if we can currently express that signature in the fully general case. If Varargs were allowed at any part of a signature and covariant then you could maybe do something like vcat(Varargs{AbstractMatrix}, SparseMatrixCSC, Varargs{AbstractMatrix}) ?

For now you could likely just list a few permutations up to a certain number of inputs? There does get to be a point where if you're concatenating N-1 dense matrices with 1 sparse matrix you may not want to promote to sparse any more, but it depends on your data sizes and there's a bit of a type stability issue.

[ViralBShah]

I thought vcat(X::Union{Matrix, SparseMatrixCSC}...) should work out ok, since the case of all dense matrix vcat will be handled by the sharper definition in the dense matrix implementation. Everything else will be sparse, so seems like it should work.

[pkofod]

You're right @ViralBShah , I must have made some sort of mistake. I had an attempt, but that seemed to cause all my hcats to be sparse. I must have a bug. The following works fine, though.

function attempt(args::Int64...)
    for arg in args
        println(arg)
    end
end

function attempt(args::Float64...)
    for arg in args
        println(arg)
    end
end

function attempt(args::Union{Int64, Float64}...)
    for arg in args
        println(Float64(arg))
    end
end

attempt(3)
attempt(3,4)
attempt(3.)
attempt(3.,4.)
attempt(3, 4.)
attempt(3., 4)