Is there a nice way to specify sums or products with excluded indices?

General Usage
1

Say I need to calculate something like

y_i = \sum_{i \neq j} \frac{1}{x_i - x_j},

Where the summation excludes identical indices. In Julia, the only way I know how to do this is with nested loops and an if statement, eg:

x = [1 2 3 4 5]
y = zeros(5)

for i in eachindex(x), j in eachindex(x)
    if i != j
        y[i] += 1/(x[i] - x[j])
    end
end

This is functional, but rather clunky to write. Is there a more concise/elegant/idiomatic way to write “hey, sum/multiply these terms over all indices, but exclude this special case”? (Note: the function above is a good example of where you might encounter this excluded index summation, since it avoids a divide-by-zero error, but that’s hardly the only application.)

2

You can use InvertedIndices.jl:

using InvertedIndices
x = [1, 2, 3, 4, 5]
x[Not(1)]                       # --> [2, 3, 4, 5]
x[1] .- x[Not(1)]               # --> [-1, -2, -3, -4]
sum(inv, x[1] .- x[Not(1)]      # --> -2.0833333

y = map(i -> sum(inv, x[i] .- x[Not(i)]), eachindex(x))

Whether it’s more elegant is up to you. I think the loop is actually fine.

3

Ok, so this didn’t appear when I searched it (figures), but did appear in the “Related Posts” section after I posted this topic. Turns out the answer was, in fact, pretty simple: How to sum every element but one - #3 by greg_plowman

4

That’s pretty slick! I might use it. Thanks for showing me!

5

To be honest, your original code was readable and made sense to me.

6

Note that this allocates a temporary array for each sum call.

An allocation-free one-liner is:

y = map(i -> sum(j -> (j != i) * inv(x[i]-x[j]), eachindex(x)), eachindex(x))

Note that this works because false is a “strong zero” in Julia:

julia> false * Inf
0.0

julia> false * NaN
0.0

Beware that the more obvious code j != i ? inv(x[i]-x[j]) : 0 is type-unstable (while using 0.0 for the other branch assumes Float64 precision rather than inferring the precision from x).

All-in-all, the “one-liner” solutions seem less clear than a nested loop, though they do have the advantage that sum uses pairwise summation for improved accuracy. map also infers the result type for you, whereas if you use zeros approach you would need to do a type computation like y = zeros(typeof(inv(oneunit(eltype(x)))), length(x)) if you want the code to be type-generic.