How much faster can we make `sum` with `skipmissing`?

[pdeffebach]

Just did some benchmarking about sum with skipmissing. It's slower than i expected, even with skipmissing's optimized reduce infrastructure.

julia> x = [rand() < .2 ? missing : rand() for i in 1:100_000];

julia> y = collect(skipmissing(x));

julia> using BenchmarkTools;

julia> @btime (sum ∘ skipmissing)($x);
  229.778 μs (6 allocations: 96 bytes)

julia> @btime (sum ∘ skipmissing)($x);
  227.911 μs (5 allocations: 80 bytes)

julia> @btime sum($y);
  12.180 μs (0 allocations: 0 bytes)

julia> @btime sum($y);
  12.185 μs (0 allocations: 0 bytes)

julia> sum2(x) = sum(skipmissing(x));

julia> @btime sum2($x);
  228.483 μs (5 allocations: 80 bytes)

julia> @btime sum2($x);
  229.343 μs (5 allocations: 80 bytes)

Is there anything to make this faster?

[KristofferC]

Is there anything to make this faster?

Is the data for the missing part guaranteed to be zero initialized? In that case you could YOLO it with:

julia> xx = unsafe_wrap(Array, Ptr{Float64}(pointer(x)), (length(x),));

julia> @btime sum($xx)
  11.961 μs (0 allocations: 0 bytes)
39871.58366123557

julia> @btime sum(skipmissing($x))
  80.376 μs (5 allocations: 80 bytes)
39871.58366123556

Probably not a good idea :P

[quinnj]

If the structure is specifically Vector{Union{T, Missing}}, where T is isbitstype, then we do guarantee the data is zero-initialized.

[vtjnash]

No, in fact, we encourage codegen to put garbage there

[quinnj]

Can you explain your comment more @vtjnash ? My comment was based on

julia/src/array.c

Line 140 in 11016bf

memset(data, 0, tot);

.

[vtjnash]

We do not guarantee it, that code is simply being inefficient. The data in memory might even be zero, but that doesn't mean the read must return zero. Reading any data from one Array type through an Array of a different type (e.g. via unsafe_wrap) is prohibited by TBAA rules, and is always permitted to introduce arbitrary undefined behavior into your results. Thus we may, at our convenience, assert in the optimizer that these pointers cannot be equal, and mis-optimize any code region that violates that assumption, or leads up to or depends upon the result of such a code region.

This appears to be a problem with the @simd macro. In particular, we mark the SkipMissing map-reduce loop, but we don't handle conditional reductions in the lowering pass, so we can't actually vectorize it. Thus with --math-mode=fast to forceable annotate it, we are much faster (2x the best instead of 20x).

[sl-solution]

To complicate matters I run another benchmark about skipmissing. ##JuliaData/DataFrames.jl#2731 (comment)

I used a "naive" approach, and I call it naive because the functions that I defined are simple in the sense that they only works for one-dimensional Floats and Ints, and the computation error is not a priority. The definition of the functions are postponed to the end of this note for clarity.

x = rand(10^6)
x = convert(Vector{Union{eltype(x),Missing}},x)
x[rand(1:length(x),1000)] .= missing

@btime sum(skipmissing($x))
  881.652 μs (6 allocations: 96 bytes)
499317.4491114656

@btime mean(skipmissing($x))
  876.515 μs (0 allocations: 0 bytes)
0.5003165813243889

@btime var(skipmissing($x))
  7.732 ms (0 allocations: 0 bytes)
0.08336800351336347

using the naive approach

@btime naive_sum($x)
  328.217 μs (0 allocations: 0 bytes)
499317.4491114656

@btime naive_mean($x)
  365.206 μs (0 allocations: 0 bytes)
0.5003165813243704

@btime naive_var($x)
  451.486 μs (0 allocations: 0 bytes)
0.08336800351336499

The function definitions

function naive_sum(f, x::AbstractArray{Union{T,Missing},1}) where {T <: Number}
    all(ismissing, x) && return missing
    _dmiss(y) = ifelse(ismissing(f(y)),  zero(T), f(y))
    mapreduce(_dmiss, +, x)
end
naive_sum(x::AbstractArray{Union{T,Missing},1}) where {T <: Number} = naive_sum(identity, x)


function naive_mean(f, x::AbstractArray{Union{T,Missing},1}) where {T <: Number}
    _op(y1,y2) = (y1[1]+y2[1],y1[2]+y2[2])::Tuple{T,Int}
    _dmiss(y) = (ismissing(f(y)) ? zero(T) : f(y), !ismissing(f(y)))::Tuple{T,Bool}
    sval, n = mapreduce(_dmiss, _op, x)::Tuple{T, Int}
    n == 0 ? missing : sval/n
end

naive_mean(x::AbstractArray{Union{T,Missing},1}) where {T <: Number} = naive_mean(identity, x)

function naive_var(f, x::AbstractArray{Union{T,Missing},1}, df=true) where {T <: Number}
    _opvar(y1,y2) = (y1[1]+y2[1], y1[2]+y2[2], y1[3]+y2[3])::Tuple{T,T,Int}
    _dmiss(y) = ismissing(f(y)) ? zero(T) : f(y)
    _varf(y) = (_dmiss(y) ^ 2, _dmiss(y) , !ismissing(f(y)))

    ss, sval, n = mapreduce(_varf, _opvar, x)::Tuple{T, T, Int}

    if n == 0
        return missing
    elseif n == 1
        return zero(T)
    else
        res = ss/n - (sval/n)*(sval/n)
        if df
            return (n * res)/(n-1)
        else
            return res
        end
    end
end

naive_var( x::AbstractArray{Union{T,Missing},1}, df=true) where {T <: Number} = naive_var(identity, x, df)

[sl-solution]

I added median to the benchmark: (increased the sample size)

x = rand(10^8)
x = convert(Vector{Union{Float64,Missing}},x)
x[rand(1:length(x),10000)] .= missing
@time median(skipmissing(x))
  4.362019 seconds (31 allocations: 912.361 MiB, 1.75% gc time)
0.5000651115865347
@time naive_median(x)
  2.353116 seconds (3 allocations: 858.307 MiB)
0.5000651115865347

[brenhinkeller]

Some updated timings, as of 1.8.3 on M1 mac

julia> versioninfo()
Julia Version 1.8.3
Commit 0434deb161e (2022-11-14 20:14 UTC)
Platform Info:
  OS: macOS (arm64-apple-darwin21.3.0)
  CPU: 10 × Apple M1 Max
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, apple-m1)
  Threads: 1 on 8 virtual cores

julia> x = [rand() < .2 ? missing : rand() for i in 1:100_000];

julia> y = collect(skipmissing(x));

julia> using BenchmarkTools

julia> @btime sum($y)
  12.791 μs (0 allocations: 0 bytes)
40151.880206376205

julia> @btime (sum ∘ skipmissing)($x)
  150.583 μs (5 allocations: 80 bytes)
40151.8802063762

Looks like some relative progress, but safe to say there's still more performance on the table.

I will also throw another option into the ring for comparison, at least for those working with AbstractFloats:

julia> z = map(i -> (ismissing(i) ? NaN : i), x);

julia> using NaNStatistics

julia> @btime nansum($z)
  19.291 μs (0 allocations: 0 bytes)
40151.88020637621

and just for fun

julia> using VectorizedStatistics

julia> @btime vsum($y)
  7.375 μs (0 allocations: 0 bytes)
40151.8802063762

so there's probably also still some more simd-ability to be had in plain base sum for that matter