Sparse: bugfix + improvements in [hv]cat + more tests

[carlobaldassi]

I fixed the bug which was affecting the sparse ref (mentioned in #585). It turned out it was actually in sparse matrix multiplication (which is used by ref).

I also fixed the promotion in sparse [hv]cat functions and added more tests.

[ViralBShah]

What is the purpose of using full() in abstractarray.jl as you have done?

[carlobaldassi]

Suppose you have a sparse matrix s and a dense matrix d.
If you do [s d] or [s; d] or [d d; s d] you get an error without full (complaining about not finding a suitable similar instantiation).
The logic here was:

• If all arguments in a [hv]cat are sparse, julia uses the sparse cat versions, which are more specialized
• If any of the arguments is not sparse, the fallback is used. And we basically want to convert everything to dense representation, hence the call to full (which doesn't do anything for dense matrices, therefore my assumption is that it would be optimized out by the compiler).

The actual drawback in performance could come when the matrix (or vector) is sparse and full is called just to be passed to similar, which is a waste. But to address this in abstractarray.jl one should define a general interface (let's call it typeof_as_full) which is guaranteed to return the outcome of typeof(full(A)) and which could be specialized for sparse arrays. Then, one could substitute similar(full(A), n, m) with Array(typeof_as_full(A), n, m) (provided that similar(A, n, m) is always equivalent to Array(typeof(A), n, m), which I'm not 100% sure about).

[carlobaldassi]

Sorry, that should have been eltype_as_full(A) set to equal eltype(full(A)). and then one could do eltype_as_full(A::SparseMatrixCSC{Tv,Ti}) = Tv

[carlobaldassi]

Maybe it is safer to add a similar_to_full (or similar_to_dense) function which is just the same as similar for non-sparse arrays, and specialize that instead.

[ViralBShah]

I think it is reasonable to convert all sparse matrices to dense in the case of sparse-dense concatenation. Perhaps we can get this to work without touching the abstractarray hvcat definitions. How about having the following in sparse.jl, which just converts all its inputs to dense and then calls the dense hvcat?

hvcat{T}(rows::(Int...), as::Union(Matrix{T},SparseMatrixCSC{T})...)

-viral

On 15-Mar-2012, at 6:39 PM, Carlo Baldassi wrote:

Suppose you have a sparse matrix s and a dense matrix d.
If you do [s d] pr [s; d] or [d d; s d] you get an error without full (complaining about not finding a suitable similar instantiation).
The logic here was:

• If all arguments in a [hv]cat are sparse, julia uses the sparse cat versions, which are more specialized
• If any of the arguments is not sparse, the fallback is used. And we basically want to convert everything to dense representation, hence the call to full (which doesn't do anything for dense matrices, therefore my assumption is that it would be optimized out by the compiler).

The actual drawback in performance could come when the matrix (or vector) is sparse and full is called just to be passed to similar, which is a waste. But to address this in abstractarray.jl one should define a general interface (let's call it typeof_as_full) which is guaranteed to return the outcome of typeof(full(A)) and which could be specialized for sparse arrays. Then, one could substitute similar(full(A), n, m) with Array(typeof_as_full(A), n, m) (provided that similar(A, n, m) is always equivalent to Array(typeof(A), n, m), which I'm not 100% sure about).

---

Reply to this email directly or view it on GitHub:
#592 (comment)

[carlobaldassi]

That's the first thing I tried indeed...
The problem in all my tests along that line was that it either keeps calling itself or the julia parser issues an ambiguity warning. That's when I gave up and patched abstractarray instead.

[ViralBShah]

Fair enough. For now, I think it is ok to just go with full(). We can do this kind of a thing whenever we have more array types. Can you just check that your full() fix does not in any way slow down the dense hvcat? I don't think it should since full() for dense arrays is essentially a no-op, but good to just be on the safe side.

-viral

On 15-Mar-2012, at 6:59 PM, Carlo Baldassi wrote:

Maybe it is safer to add a similar_to_full (or similar_to_dense) function which is just the same as similar for non-sparse arrays, and specialize that instead.

---

Reply to this email directly or view it on GitHub:
#592 (comment)

[carlobaldassi]

I just did some tests and found no measurable difference. This is the code I used:

println("hcat")
tic()
for i = 1 : 100000
    a = rand(10, 10)
    b = [a a]
    c = [a a' a]
end
toc()

println("vcat")
tic()
for i = 1 : 100000
    a = rand(10, 10)
    b = [a; a]
    c = [a; a'; a]
end
toc()

println("hvcat")
tic()
for i = 1 : 100000
    a = rand(10, 10)
    b = [a a; a a]
    c = [a a' a; a' a a; a a a']
end
toc()

[carlobaldassi]

I did some more tests with more samples and smaller matrices (to emphasize the function call overhead) and this time fixing the random seed, and again I don't see any measurable difference in performance (i.e. less than statistical fluctuations across many runs).
I'd like to do the same with sparse matrices and a tentative implementation of similar_as_full now.

[carlobaldassi]

I did the tests implementing similar_to_full. There doesn't seem to be any condition of matrix size and density where a measurable performance gain can be achieved (I tried up to 1000x1000 matrices). The only benefit would come from sparing some temporary memory allocation; I'm not sure if this is worth the burden of carrying over another interface for AbstractArrays.

[ViralBShah]

I suggest let's just leave it the way it is for now. The current code is more readable, and we are not going to require all this flexibility any time soon. Let's focus our effort on adding more sparse matrix functionality. :-)

-viral

On 15-Mar-2012, at 10:33 PM, Carlo Baldassi wrote:

I did the tests implementing similar_to_full. There doesn't seem to be any condition of matrix size and density where a measurable performance gain can be achieved (I tried up to 1000x1000 matrices). The only benefit would come from sparing some temporary memory allocation; I'm not sure if this is worth the burden of carrying over another interface for AbstractArrays.

---

Reply to this email directly or view it on GitHub:
#592 (comment)