RFC: Add weights argument to sum

[nalimilan]

This will be consistent with the weights argument to be added to similar functions in Statistics. Since that code is needed notably to computed the weighted mean in Statistics, I figured it is better to make it public rather than only defining internal helpers for mean.

This is based on code from StatsBase, of which one part needs to live in LinearAlgebra since it depends on BLAS and dot.

Extracted from #31395.

[stevengj]

Why should weights have to be an array, as opposed to just an iterable?

[nalimilan]

That's how things work currently in StatsBase. I guess we could make this more general at least for the simple case where dims=:. But that would require adding a separate method as @simd only works when using indices.

[stevengj]

I wonder whether this should be an iterator instead? i.e. you could imagine:

sum(Iterators.multiply(x, weights))

[stevengj]

This will be much less accurate than the non-weighted algorithm based on pairwise summation.

[nalimilan]

Yes, probably (I haven't considered this, that's just a port of the StatsBase code with some cleanup). We could improve this, but note that for BlasReal the BLAS-optimized dot is called instead, which has different properties anyway (not sure whether it's going to be more accurate or less).

If we want to preserve accuracy, that would be an argument in favor of this PR, since a naive broadcast call won't use pairwise summation.

[stevengj]

In particular, I could imagine an algorithm in which Iterators.multiply(x, weights) produces either:

1. An Iterators.Multiply(x, weights) object, similar to z[1]*z[2] for z in zip(x,weights) except having an eltype (allowing e.g. empty sums to work correctly).

2. If x and weights are both abstract arrays, then a MultipliedArray(x, weights) <: AbstractArray object for which getindex(_, i) produces x[i]*weights[i].

Both of these could then be passed to the existing sum function (as well as other functions, e.g. external packages like xsum), and in the second case it would do an accurate pairwise summation.

(It could also be made even more generic by passing a binary op(x,y) to the iterator.)

Seems like it would take less code than what you are doing now, would be more accurate (for arrays), and would be more composable.

[stevengj]

Or we could just wait for sum(@lazy x .* weights) ala #19198.

[stevengj]

But with regards to adding this functionality to sum specifically, why do we need that when we have dot(weights, x) or @views dot(weights, x[...]), since it seems like that is basically the functionality here?

[simeonschaub]

@stevengj Wouldn't sum(Base.broadcasted(*, x, weights)) do the exact same thing?

[nalimilan]

As I wrote in the description, the main point of this PR is that we're going to add the weights keyword argument to many statistics functions (mean, var, cor, quantile...), so adding it to sum would be consistent. I'm not too attached to it, we can keep it internal for use by mean in Statistics.

Another reason is that the current behavior (implemented in StatsBase) differs from broadcast in that you provide a weights vector corresponding to each slice defined according to the dims argument: row weights when summing over columns, or vice-versa. So the Iterators.multiply(x, weights) or @lazy x .* weights approach won't be enough to replicate this behavior. I can't say how useful it is for sum (and mean), but that sounds logical, and again consistent with what is needed for e.g. cor (where an observation is a slice, not a single value).

A final reason can be that the API from this PR allows using efficient and/or accurate algorithms under the hood, e.g. pairwise summation or BLAS-powered dot (as you noted above).

If we decide that the weighted sum/mean over dimensions use case isn't worth all that code, then yes, this PR isn't very useful.

[tkf]

pairwise summation

This would be handled if we can use mapreduce on Broadcasted, right? This is implemented in #31020. It would be nice if we can reuse the code in #31020 here even if it were just for internals (#31020 is not about the surface syntax).

BLAS-powered dot

In principle, we can handle sum(@lazy x .* y) specially to call BLAS.dot directly if x and y are both DenseArray etc., because the information required to do this is accessible in the type of the Broadcasted object. Not sure if it is worth the trouble, tough.

[nalimilan]

This [pairwise summation] would be handled if we can use mapreduce on Broadcasted, right?

Ah, yeah, I guess we could add special methods for that. Not sure whether that's OK given that it would give different results from the general algorithm, but since increased accuracy is a good thing maybe that's not a problem. The same applies to using dot.

[oxinabox]

We should have this.
It makes sense especially as sum in the generalization of count

If one wants a weighted tollerance for a cerain amount of missing data then one can right now get weighted relative confidance like:

mean(ismissing, xs, weights) < 0.1|| error("too much missing data > 10% under weighting")

But if you want to do this with an absolute count,
say 10 weighted missing values,
you can't.

sum would fix that.

---

This generalizes to over dimentions.
Real code I am using to drop rows will too much missing

    missings = ismissing.(X)  # no `mean(f, x, wv; dims)` function
    mask = vec(mean(missings, wv..., dims=1) .< 0.1)
    X[mask, :]

Currently harder to write the absolute tolerance version

[StefanKarpinski]

My worry is that this opens the door to adding weights support everywhere. Why can't this be done externally in say a Weights.jl package? The methods seem like they should be non-overlapping.

[nalimilan]

Well keyword arguments cannot be handled in packages unfortunately...

The current design in StatsBase relies on dispatch, with methods like sum(::AbstractArray, ::AbstractWeights), requiring people to do things like sum(x, weights(w)). That gives some weird things sometimes, e.g. var(x, weights(w)) is the weighted variance of x but cov(x, weights(w)) is the covariance of x and w, since w isa AbstractWeights <: AbstractArray (JuliaStats/StatsBase.jl#409).

Also, people have often complained that they expected e.g. sample(x, w) to work, but instead they had to write sample(x, weights(w)) (JuliaStats/StatsBase.jl#335). A keyword argument would be more standard and probably less confusing. Overall the problem is that for some signatures there's no ambiguity so we could accept any AbstractVector as weights, but for others we have to require AbstractWeights for dispatch. In the end the result is a big inconsistent and confusing. Probably not the end of the world though. @ararslan will likely want to comment as he was a big supporter of the keyword argument approach.

[oxinabox]

I also support the kwarg approach.
It just feels right, weights are like dims.

[tkf]

but cov(x, weights(w)) is the covariance of x and w

Maybe it was already considered, but why not use cov(weighted(x, w)) where weighted(x, w) is effectively a lazy x .* w?

[StefanKarpinski]

weighted(x, w) seems much more composable and can be implemented completely externally.

[vtjnash]

theoretically, sum(for x .* w) is even more composable (and concise)

[StefanKarpinski]

I disagree: it's often convenient to associate values and weights together a single time and thereafter use them as though they were just values (implicitly weighted). Using a lazy construct like for x .* w can be consumed only once and doesn't support indexing, etc.

[nalimilan]

As noted before, broadcasting wouldn't work as the weights vector is broadcast to the dimension over which reduction in performed, rather than always treated as a column vector.

But weighted(w, x) is interesting. However, note that for operations on multiple arguments it would differ from the unweighted syntax: cov(x, y) would become cov(weighted(w, x, y)). I guess that would be OK if we consider that weighted does something like hcat. But if we really want it to return an AbstractArray object it will have to take an additional argument specifying whether weights apply to rows or columns when the argument is a matrix, which will imply a repetition like sum(weighted(w, x, dims=1), dims=1), unless we want to automatically assume that you want to reduce along that dimension (which is the most useful operation). Need to think about it.

[pdeffebach]

Is this getting merged? I support using keyword arguments instead of types, in which case this would need to live here as opposed to Statistics.

[nalimilan]

For the record, something I hadn't realized when I wrote my comments above: weighted(w, x) absolutely cannot be a lazy equivalent of w .* x, as weighting is equivalent to multiplication by weights only in very peculiar cases like sum and mean. For quantile you certainly don't want to multiply each entry by its weight, but rather repeat it (at least for frequency weights). For var, the bias correction depends on the kind of weight. And so on.

So the wrapper approach would only be used for dispatch, and the resulting object would implement almost no methods. Each function would have to opt-in to support it as appropriate in its particular case.

[ViralBShah]

Given that we want to move Statistics out (#46501), perhaps we should close this and adopt what we need in Statistics.jl?