Error in mean value of sparse arrays

in the following code, I noticed two errors:

  1. the mean value of the positive elements is less than the mean value including zeros
  2. the mean value obtained from the loop is different than the mean value from the function “mean”.
n = 5000
sp_mat = sprand(Float32,n,n,0.99)./20
for i=1:n
	for j=1:n

I added some more details showing what is computed:

using SparseArrays
using Statistics

n = 5000
sp_mat = sprand(Float32,n,n,0.99)./20
@show mean(sp_mat)
@show sum(nonzeros(sp_mat)) / (size(sp_mat, 1) * size(sp_mat, 2))
@show mean(sp_mat[sp_mat.>0])
@show size(sp_mat[sp_mat.>0]), sum(sp_mat[sp_mat.>0]) / size(sp_mat[sp_mat.>0], 1)
for i=1:n
	for j=1:n
		global sum_loop+=sp_mat[i,j]
@show sum_loop/n/n


mean(sp_mat) = 0.02475363f0
sum(nonzeros(sp_mat)) / (size(sp_mat, 1) * size(sp_mat, 2)) = 0.02475363f0
mean(sp_mat[sp_mat .> 0]) = 0.025003498f0
(size(sp_mat[sp_mat .> 0]), sum(sp_mat[sp_mat .> 0]) / size(sp_mat[sp_mat .> 0], 1)) = ((24750166,), 0.025003498f0)
(sum_loop / n) / n = 0.023831278f0

I believe the reason your loop is showing a different result is due to sum using pairwise summation.

1 Like

And of course, this one:


is wrong as the number of elements for mean is smaller than in sp_mat.

1 Like

This is the last routine which does the sum:

@noinline function mapreduce_impl(f, op, A::AbstractArrayOrBroadcasted,
                                  ifirst::Integer, ilast::Integer, blksize::Int)
    if ifirst == ilast
        @inbounds a1 = A[ifirst]
        return mapreduce_first(f, op, a1)
    elseif ifirst + blksize > ilast
        # sequential portion
        @inbounds a1 = A[ifirst]
        @inbounds a2 = A[ifirst+1]
        v = op(f(a1), f(a2))
        @simd for i = ifirst + 2 : ilast
            @inbounds ai = A[i]
            v = op(v, f(ai))
        return v
        # pairwise portion
        imid = (ifirst + ilast) >> 1
        v1 = mapreduce_impl(f, op, A, ifirst, imid, blksize)
        v2 = mapreduce_impl(f, op, A, imid+1, ilast, blksize)
        return op(v1, v2)

In case of SparseArrays the above function is called with a view over only the non-zero values, which is SparseArrays.nzvalview(sp_mat).

You can see, that there is no difference now, by calling it directly:

julia> using SparseArrays, Statistics

julia> n = 5000

julia> sp_mat = sprand(Float32,n,n,0.99)./20;

julia> sum(sp_mat)

julia> inds = LinearIndices(SparseArrays.nzvalview(sp_mat));

julia> Base.mapreduce_impl(identity, +, SparseArrays.nzvalview(sp_mat), first(inds), last(inds), Base.pairwise_blocksize(identity, +))

The sum algorithms sums up blockwise, with blocksize of 1024. (better short description below)
Summing up floating point numbers in different orders gives different results because of floating point precision.
You can find out which function is doing the actual calculation by using @which and @less with the next function you see in @less, e.g.:

julia> @less sum(sp_mat)
@inline ($fname)(a::AbstractArray; dims=:, kw...) = ($_fname)(a, dims; kw...)

julia> @less Base._sum(sp_mat,:)
($_fname)(a, ::Colon; kw...) = ($_fname)(identity, a, :; kw...)
julia> @less Base._sum(identity,sp_mat,:)
($_fname)(f, a, ::Colon; kw...) = mapreduce(f, $op, a; kw...)
julia> @which Base._sum(identity,sp_mat,:) 
#because we need to see how $op is defined, which is :add_sum
_sum(f, a, ::Colon; kw...) in Base at reducedim.jl:878

julia> @less mapreduce(identity,Base.add_sum,sp_mat,:)
mapreduce(f, op, itrs...; kw...) = reduce(op, Generator(f, itrs...); kw...)

And so on until you reach the worker function which does the actual work.
I like to do this sometimes, despite it can be tedious, because there is always so much to learn.

1 Like

More precisely, the default sum adds pairwise recursively, with a base case of 1024 where it switches to a simple loop.


thanks, this is amazingly informative. If I run the following, I get

  1. different values using “sum” and the loop
  2. different values if I run the loop for different permutations.

So yes, I guess it was a combination of a different order of summing up, and of course Float32 precision along with small values to be added.

using Random
rand_vec = rand(Float32,nn)./20
sum_loop = Float32(0)
ii = randperm(nn)
for i in ii
   sum_loop += rand_vec[i]