Fastest quadratic mixing rule?

General Usage Performance
1

I’m in need of a good mixing rule function. A mixing rule is a type of of weighed mean. if we have a property p and some weights x, a mixing rule is like this:

p_mixing = sum(op(pi,pj)*xi*xj for i,j = 1:N)

the op(pi,pj) is a type of mean operation, where op(pi,pi)= pi (for example, an aritmethic mean)

using the properties of symmetry and that fact, i made a function that does such mixing rule

function mixing_rule(op,x,p)
    N = length(x)
    res1 = p[N]*x[N]^2
    for i = 1:N-1
        res1 += p[i]*x[i]^2
        for j=(i+1):N
            res1 += 2*x[i]*x[j]*op(p[i],p[j])
        end
    end
    return res1
end

can it be faster?, i’m not convinced of my methods, or, there exists some function in some package that does this? sounds like it can be implemented by a mapreduce or somethig

2

This takes less than 1/6th the time on my machine:

function mixing_rule_2(op,x,p)
    N = length(x)
    @boundscheck length(p) == N || throw(DimensionMismatch("Dimensions of x and p don't match"))
    @inbounds begin
        res1 = zero(p[1] * x[1]^2)
        for i = 1 : N
            res1 += p[i] * x[i]^2
            @simd for j = 1 : i - 1
                res1 += 2 * x[i] * x[j] * op(p[i], p[j])
            end
        end
    end
    return res1
end

  • use @inbounds and @simd to elide bounds checks and get the inner loop to vectorize

  • iterate over j = 1 : i - 1 instead of j=(i+1):N in the inner loop for better cache locality

You can probably do quite a bit better still (e.g. with threading), but this is a nice quick win.

Full code including benchmark setup. 
function mixing_rule(op,x,p)
    N = length(x)
    res1 = p[N]*x[N]^2
    for i = 1:N-1
        res1 += p[i]*x[i]^2
        for j=(i+1):N
            res1 += 2*x[i]*x[j]*op(p[i],p[j])
        end
    end
    return res1
end

function mixing_rule_2(op,x,p)
    N = length(x)
    @boundscheck length(p) == N || throw(DimensionMismatch("Dimensions of x and p don't match"))
    @inbounds begin
        res1 = zero(p[1] * x[1]^2)
        for i = 1 : N
            res1 += p[i] * x[i]^2
            @simd for j = 1 : i - 1
                res1 += 2 * x[i] * x[j] * op(p[i], p[j])
            end
        end
    end
    return res1
end

using BenchmarkTools
N = 10_000
x = rand(N)
p = rand(N)
op = (x, y) -> (x + y) / 2

@assert mixing_rule(op, x, p) ≈ mixing_rule_2(op, x, p) atol=1e-5
@btime mixing_rule_2($op, $x, $p)

@btime mixing_rule($op, $x, $p)
4 Likes
3

Thanks!, it worked flawlessly, to optimize further, i found this:

@boundscheck checkbounds(p,N)

and it shaves some aditional time, but nothing compared to your original function :grin: