# Optimizing sums of products (dot products)

**URL:** https://discourse.julialang.org/t/optimizing-sums-of-products-dot-products/119515
**Category:** Performance
**Tags:** simd, linearalgebra, sum
**Created:** [September 17, 2024, 4:40pm UTC](https://discourse.julialang.org/t/optimizing-sums-of-products-dot-products/119515 "2024-09-17T16:40:15Z")
**Posts on this page:** 1
**Showing post:** 11

<div class="post-metadata">

### Author: ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)
#### Post date: [September 23, 2024, 5:50pm UTC](https://discourse.julialang.org/t/optimizing-sums-of-products-dot-products/119515/11 "2024-09-23T17:50:09Z")

</div>

I’m always on the hunt for good “mere mortal” code that goes awry.

One of the fundamentals is that fast math _allows the compiler to **sometimes** give you different answers for basic maths_ depending on the context and surrounding code and available optimizations. That makes debugging very challenging, and especially so in Julia where small functions often inline and constant propagate and then subsequently optimize in different ways.

So a huge practical consequence of this — yes even in “mere mortal” higher-level less-numericsy code — is that you can’t assume that the maths you did in an `if` conditional will match the results of the **same maths** elsewhere. Just because `f() < g()` in one place in your code doesn’t necessarily mean that `f() != g()` with ffast/funsafe-math. In fact, `f()` might be bigger than `g()` elsewhere. So even if you try to “defend” against fast-math with branches, your adversary is _already there in your defense_, mucking around in your conditionals. That most infamously happens with `isnan` and `isinf` checks (and @mikmoore has some great examples above), but it also applies to finite math, too — subtractions can get re-arranged to be additions, things can swap sides of the equality, etc., etc.

My favorite examples are in the `contract` category because it feels like converting a mul-add to the **higher-precision** fused multiply and add should be a no-brainer — better accuracy _and_ better performance, what’s not to love?

- The expression `a*b + c*d` now has _three possible answers_, and the one that’s most accurate is value-dependent. This is particularly bad for complex multiplication, because `x*conj(x)` definitely needs to be real — and it is, but only if you **don’t contract**. It’s actually quite relevant to the original conversation here, because that’s precisely what’s on the inside of the dot-product `x'x`. With FMAs, the imaginary parts would be sometimes positive, sometimes negative. Whatever, Matt, you might say, I’m not implementing complex arithmetic…
- Taking the square root of a negative float is an error, so you might have a check to ensure that your argument is positive in that branch. What if you’re solving the quadratic equation? The operand to `sqrt` there is `b*b - 4*a*c` — look familiar? Just like `a*b + c*d`, there are multiple possible ways to compute this with FMAs in the mix. I’d love to find time to build a realistic example that demonstrates different values in an `if` than the `sqrt`, but it’s the same idea as [this example](https://discourse.julialang.org/t/when-if-a-b-x-1-a-b-divides-by-zero/7154/5). What would make it even more fun is that — upon getting a spurious domain error — you might add logging and then _never see it again_, because just looking at the intermediate values would break the `fma`.

---

_[View the full topic](https://discourse.julialang.org/t/optimizing-sums-of-products-dot-products/119515)._