Wrong weighted sum and mean results for complex numbers

[aplavin]

Weighted mean gives wrong (conjugated) result for complex arrays: mean([1 + 0.5im], weights([1])) == 1.0 - 0.5im.

[nalimilan]

Indeed. That's because we use dot([1 + 0.5im], [1]) to compute the weighted sum (for performance reasons). I guess we could simply do dot([1], [1 + 0.5im]) since only the first argument is conjugated, and we know weights are real.

Would you make a pull request to apply that change (in weights.jl) and add a test?

[oschulz]

Just ran into this in another context: #534 isn't quite the same, because it's about values that are vectors of real numbers, but the cause is the (recursive) behavior of dot(...) as well.

I think using dot in sum(A, w) is wrong, semantically. Could we change the definition of sum(A, w) to something like sum(A .* values(w)) or mapreduce(*, +, A, values(w))? I can do a PR, if there are no objections.

[nalimilan]

The reason we use dot AFAICT is that it's much more efficient thanks to BLAS. It's not an issue to abuse it if we are sure it's correct in all cases in which we use it. What we can do is to have a method using it for arrays of numbers, and a fall back method using a slower but generally correct approach.

[oschulz]

Sounds good. Regarding on where do draw the line - BLAS would only accelerate Real, right (resp. do the wrong thing for Complex for this use case)?

So we'd do

Base.sum(v::AbstractArray{<:Real}, w::AbstractWeights) = dot(v, values(w))

Base.sum(v::AbstractArray, w::AbstractWeights) = ...

right? @nalimilan, what do you think would be best for the generic implementation:

# Option 1:
sum(A .* values(w))

# Option 2:
sum(Base.Broadcast.Broadcasted(*, (A, values(w))))

# Option 3:
mapreduce(*, +, A, values(w))

[oschulz]

I've been thinking - maybe we should drop use of dot altogether, though, even if it means a performance loss, since we'd gain precision. AFAIK, Julia's sum() uses pairwise summation, to keep sums stable when summing over values with large differences in magnitude. But I don't think BLAS does fancy summing:

julia> A = Float32.(vcat(1.0, fill(1E-10, 10^6)))
julia> w = Weights(fill(Float32(1.0), length(A)))

julia> sum(A)
1.0000999f0

julia> sum(A, w), dot(A, values(w))
(1.0000969f0, 1.0000969f0)

julia> mapreduce(*, +, A, values(w))
1.0000999f0

Though I just saw that we're using a simple sum in StatsBase._moment2 instead of Welford's algorithm or similar - is that intentional? Statistics._varm uses centralize_sumabs2 for stability.

[nalimilan]

Options 2 and 3 are the best ones. Not sure how to choose between them, I guess benchmarking is needed.

mapreduce(*, +, A, values(w)) doesn't use pairwise summation either, so I think using BLAS is better. Can you check whether sum(Base.Broadcast.Broadcasted(*, (A, values(w)))) use pairwise summation?

IIRC we use Welford's algorithm when possible in var, but not always. Please file a separate issue if that's not the case.

[oschulz]

mapreduce(*, +, A, values(w)) doesn't use pairwise summation either [...]

Yes, it does (AFAIK the default Base.sum is actually defined via `mapreduce, which operates pairwise). Here's a more extreme example with a significant differenct in result:

julia> A = vcat([Float32(1E0)], fill(Float32(1E-8), 10^8));

julia> w = Weights(fill(Float32(1.0), length(A)));

julia> sum(A)
1.9999989f0

julia> mapreduce(*, +, A, w)
1.9999989f0

julia> sum(A, w)
1.9384409f0

julia> dot(A, values(w))
1.9384409f0

IIRC we use Welford's algorithm when possible in var, but not always. Please file a separate issue if that's not the case.

Oh, sure - I just stumbled over it while looking into sum and wasn't sure if it was intentional.

[nalimilan]

Actually mapreduce(*, +, A, w) allocates a copy to store the result, so that's not an option.

[oschulz]

Actually mapreduce(*, +, A, w) allocates a copy to store the result, so that's not an option.

Oops, indeed - I naively thought it would be allocation-free. In that case, mapreduce is out, of course.

[oschulz]

sum(Base.Broadcast.Broadcasted(*, (A, values(w)))) seems to be allocation-free, but it's fairly slow.

[oschulz]

Ok, so maybe we do

Base.sum(v::AbstractArray{<:Real}, w::AbstractWeights) = dot(v, values(w))

Base.sum(v::AbstractArray, w::AbstractWeights) = sum(Base.Broadcast.Broadcasted(*, (A, values(w))))

for now, to get correct behavior in the general case without performance loss in the specific case (real numbers), and worry about precision improvement in a separate issue? Would that be fine with you, @nalimilan?

[nalimilan]

As I noted in my first comment, I think this will work:

Base.sum(v::AbstractArray{<:Complex}, w::AbstractWeights) = dot(w, v)

For the general case, it looks like broadcasting is slower than the dot fallback:

julia> x = rand(10_000);

julia> z = Vector{Union{Float64, Missing}}(x);

julia> @btime dot(x, z);
  13.443 μs (1 allocation: 16 bytes)

julia> f(A, w) = sum(Base.Broadcast.Broadcasted(*, (A, w)))
f (generic function with 1 method)

julia> @btime f(x, z);
  22.258 μs (11 allocations: 240 bytes)

Since we don't actually need the broadcasting features, I guess we should simply use a loop similar to dot, but without calling dot recursively:

function dot(x::AbstractArray, y::AbstractArray)
    lx = length(x)
    if lx != length(y)
        throw(DimensionMismatch("first array has length $(lx) which does not match the length of the second, $(length(y))."))
    end
    if lx == 0
        return dot(zero(eltype(x)), zero(eltype(y)))
    end
    s = zero(dot(first(x), first(y)))
    for (Ix, Iy) in zip(eachindex(x), eachindex(y))
        @inbounds s += dot(x[Ix], y[Iy])
    end
    s
end