# Fastest quadratic mixing rule?

**URL:** <https://discourse.julialang.org/t/fastest-quadratic-mixing-rule/22785>\
**Category:** Performance\
**Created:** [April 5, 2019, 8:32am UTC](https://discourse.julialang.org/t/fastest-quadratic-mixing-rule/22785 "2019-04-05T08:32:04Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![longemen3000](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/longemen3000/32/7298_2.png) [@longemen3000](https://discourse.julialang.org/u/longemen3000)\
**Post date:** [April 5, 2019, 8:32am UTC](https://discourse.julialang.org/t/fastest-quadratic-mixing-rule/22785/1 "2019-04-05T08:32:04Z")

</div>

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

```julia
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

---

<div class="post-metadata">

**Author:** ![tkoolen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkoolen/32/1603_2.png) [@tkoolen](https://discourse.julialang.org/u/tkoolen)\
**Post date:** [April 5, 2019, 9:48am UTC](https://discourse.julialang.org/t/fastest-quadratic-mixing-rule/22785/2 "2019-04-05T09:48:30Z")

</div>

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

```julia
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.**
>
> ```julia
> 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)
> 
> ```

---

<div class="post-metadata">

**Author:** ![longemen3000](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/longemen3000/32/7298_2.png) [@longemen3000](https://discourse.julialang.org/u/longemen3000)\
**Post date:** [April 6, 2019, 12:55am UTC](https://discourse.julialang.org/t/fastest-quadratic-mixing-rule/22785/3 "2019-04-06T00:55:07Z")

</div>

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

```julia
@boundscheck checkbounds(p,N)

```

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