Map for sparse matrices

[jumutc]

Please consider the code below where julia is converting sparse matrix into dense one and doing map on the top. Should we have map defined only for non-zero elements in sparse matrix instead? Otherwise the slow down is an issue and one cannot apply effectively λ-calculus in this case!

julia> test = sprand(10^3,10^3,.01)
1000x1000 sparse matrix with 10144 Float64 entries:
    [80  ,    1]  =  0.993039
    [128 ,    1]  =  0.117601
    [152 ,    1]  =  0.974119
    [259 ,    1]  =  0.0362442
    [289 ,    1]  =  0.621536
    [371 ,    1]  =  0.653076
    [631 ,    1]  =  0.131718
    ⋮
    [439 , 1000]  =  0.0621062
    [538 , 1000]  =  0.109039
    [613 , 1000]  =  0.212955
    [620 , 1000]  =  0.147798
    [640 , 1000]  =  0.479203
    [702 , 1000]  =  0.88309
    [884 , 1000]  =  0.780324
    [892 , 1000]  =  0.0164652

julia> f = v -> v+1
(anonymous function)

julia> @time map(f,test)
elapsed time: 14.162436633 seconds (89388704 bytes allocated)
1000x1000 sparse matrix with 1000000 Float64 entries:
    [1   ,    1]  =  1.0
    [2   ,    1]  =  1.0
    [3   ,    1]  =  1.0
    [4   ,    1]  =  1.0
    [5   ,    1]  =  1.0
    [6   ,    1]  =  1.0
    [7   ,    1]  =  1.0
    ⋮
    [993 , 1000]  =  1.0
    [994 , 1000]  =  1.0
    [995 , 1000]  =  1.0
    [996 , 1000]  =  1.0
    [997 , 1000]  =  1.0
    [998 , 1000]  =  1.0
    [999 , 1000]  =  1.0
    [1000, 1000]  =  1.0

julia> @time (I,J,V)=findnz(test); sparse(I,J,map(f,V))
elapsed time: 0.000144994 seconds (243784 bytes allocated)
1000x1000 sparse matrix with 10144 Float64 entries:
    [80  ,    1]  =  1.99304
    [128 ,    1]  =  1.1176
    [152 ,    1]  =  1.97412
    [259 ,    1]  =  1.03624
    [289 ,    1]  =  1.62154
    [371 ,    1]  =  1.65308
    [631 ,    1]  =  1.13172
    ⋮
    [439 , 1000]  =  1.06211
    [538 , 1000]  =  1.10904
    [613 , 1000]  =  1.21296
    [620 , 1000]  =  1.1478
    [640 , 1000]  =  1.4792
    [702 , 1000]  =  1.88309
    [884 , 1000]  =  1.78032
    [892 , 1000]  =  1.01647

[mlubin]

Defining map to only work on the stored elements (note: this may include zeros) of a sparse matrix seems a bit too magical and liable to lead to surprising behavior. I think it's better to be explicit in this case. You could do map!(f,nonzeros(M)) to update the stored elements in place.

[jumutc]

Maybe one could imagine something like mapnz to make it work? Also I thing this applies to other methods like filter, fold!

[JeffBezanson]

Fully agree with @mlubin.
filter is easier since it doesn't change any values, only selects from them. However it doesn't make much sense for 2-d arrays.

[StefanKarpinski]

If we had a nmatrix type with non-zero default, this would be straightforward. But we don't and no one seems that interested in it. As it stands, we could make it an error to map a function to a sparse array that doesn't map zero to zero. In that case, we could require specifying a dense result array type.

[ViralBShah]

I agree with @mlubin 's solution too. Otherwise, it is too magical for anyone else reading the code.

[JeffBezanson]

@StefanKarpinski my vague sense is that zero is special, in that common operations like matrix multiply and scaling generally (though not always) preserve it. The next generalization is to different semirings, which I believe we support via the sparse code's use of zero(T).

[jumutc]

For sure there is no place for magical transformations and behaviour. @StefanKarpinski might be right in throwing an error if one maps zeros to something else on SparseMatrixCSC object. From my side it would be nice to have some straightforward functionality like mapnz or foldnz which would work solely with nonzeros.

[JeffBezanson]

Ah, I see we also use x == 0. I guess this is kind of ambiguous if x is from an alternate semiring on the reals. I'd hate to have to use x == zero(x) everywhere though.

[nalimilan]

Raising an error would be painful as it would mean functions working with any AbstractArray would not work with sparse arrays. If you're OK to pay the (memory) price of getting a dense result, why prevent you from doing that? Else, if you're forced to convert the sparse array to a dense one first, you'll need twice as much memory (original + result).

[jumutc]

Maybe redefining the behaviour of the default map is not a good idea but extending the Base with some additional functionality which would help to work with sparse matrices seems to me reasonable.

[mlubin]

I'm not sure if a map function applied to only stored elements is even meaningful in general given that zero elements themselves could be stored. The result could change depending on stored zeros, which isn't supposed to happen. I think something like this should first live in a package before we consider adopting it in Base.

[jumutc]

But storing zeroes in sparse matrix format is jet another issue. A[i,j] = 0 for sparse matrices shouldn't lead to any modification of nzval field or other fields of SparseMatrixCSC object.

[StefanKarpinski]

@StefanKarpinski my vague sense is that zero is special, in that common operations like matrix multiply and scaling generally (though not always) preserve it. The next generalization is to different semirings, which I believe we support via the sparse code's use of zero(T).

I've outlined before how one can pretty easily have non-zero default sparse matrices decomposed as S = Z + c where Z is a zero-default sparse matrix and c is a constant. Then you can do many computations pretty easily, e.g. (Z1 + c1) + (Z2 + c2) = (Z1 + Z2) + (c1 + c2) and (Z1 + c1)*(Z2 + c2) = Z1*Z2 + c1*Z2 + c2*Z1 + c1 + c2. It's mildly complicated, but really not all that bad. The advantage is that this kind of sparse matrix is closed under all kinds of operations, including map, whereas zero-default sparse matrices are not.

[JeffBezanson]

Ah yes, I remember that now.

[kmsquire]

On Wed, May 28, 2014 at 9:10 AM, Stefan Karpinski
notifications@github.comwrote:

@StefanKarpinski https://github.com/StefanKarpinski my vague sense is
that zero is special, in that common operations like matrix multiply and
scaling generally (though not always) preserve it. The next generalization
is to different semirings, which I believe we support via the sparse code's
use of zero(T).

I've outlined before how one can pretty easily have non-zero default
sparse matrices decomposed as S = Z + c where Z is a zero-default sparse
matrix and c is a constant. Then you can do many computations pretty
easily, e.g. (Z1 + c1) + (Z2 + c2) = (Z1 + Z2) + (c1 + c2) and (Z1 +
c1)_(Z2 + c2) = Z1_Z2 + c1_Z2 + c2_Z1 + c1 + c2. It's mildly complicated,
but really not all that bad. The advantage is that this kind of sparse
matrix is closed under all kinds of operations, including map, whereas
zero-default sparse matrices are not.

I don't remember this, but I like it!

Kevin