# Fastest algorithm for (xi - xj) \* (yi - yj) for i in 1:n, j in 1:n

**URL:** <https://discourse.julialang.org/t/fastest-algorithm-for-xi-xj-yi-yj-for-i-in-1-n-j-in-1-n/106988>\
**Category:** Performance\
**Tags:** question, performance\
**Created:** [December 1, 2023, 11:19am UTC](https://discourse.julialang.org/t/fastest-algorithm-for-xi-xj-yi-yj-for-i-in-1-n-j-in-1-n/106988 "2023-12-01T11:19:14Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [December 1, 2023, 11:19am UTC](https://discourse.julialang.org/t/fastest-algorithm-for-xi-xj-yi-yj-for-i-in-1-n-j-in-1-n/106988/1 "2023-12-01T11:19:14Z")

</div>

Suppose we are given two “big” vectors `x` and `y` of real numbers (e.g. `Float64`), both with the same length `n`. What is the fastest algorithm to compute

```julia
[(x[i] - x[j]) * (y[i] - y[j]) for i in 1:n, j in 1:n]

```

?

Here is a baseline:

```julia
using BenchmarkTools

n = 10_000
x = rand(n)
y = rand(n)

function baseline(x, y)
  n = length(x)
  D = Matrix{Float64}(undef, n, n)
  @inbounds for j in 1:n
    xj = x[j]
    yj = y[j]
    for i in 1:n
      xi = x[i]
      yi = y[i]
      D[i,j] = (xi - xj) * (yi - yj)
    end
  end
end

@btime baseline(x, y)
  168.321 ms (2 allocations: 762.94 MiB)

```

Appreciate any help with the following questions:

1. Can you find a faster algorithm that exploits parallel hardware (e.g. GPU)?
2. Is it possible to rely on FFT algorithms to compute the output efficiently?
3. Given the proposed algorithm, what is the largest `n` assuming 32Gb of RAM?

---

<div class="post-metadata">

**Author:** ![heliosdrm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/heliosdrm/32/3851_2.png) [@heliosdrm](https://discourse.julialang.org/u/heliosdrm)\
**Post date:** [December 1, 2023, 11:32am UTC](https://discourse.julialang.org/t/fastest-algorithm-for-xi-xj-yi-yj-for-i-in-1-n-j-in-1-n/106988/2 "2023-12-01T11:32:40Z")

</div>

A very simple observation, although it does not go into the “difficult” parts of the question: since the resulting matrix is symmetric, you can save half of the time and memory by computing only the upper or lower triangle.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [December 1, 2023, 11:33am UTC](https://discourse.julialang.org/t/fastest-algorithm-for-xi-xj-yi-yj-for-i-in-1-n-j-in-1-n/106988/3 "2023-12-01T11:33:35Z")

</div>

Yes, I forgot to add this optimization to the baseline. Doing it now, thanks.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [December 1, 2023, 11:38am UTC](https://discourse.julialang.org/t/fastest-algorithm-for-xi-xj-yi-yj-for-i-in-1-n-j-in-1-n/106988/4 "2023-12-01T11:38:48Z")

</div>

Actually, my attempt to exploit symmetry in row-major orientation led to slower run time. Appreciate if someone can help.

---

<div class="post-metadata">

**Author:** ![Tarny\_GG\_Channie](https://avatars.discourse-cdn.com/v4/letter/t/3bc359/32.png) [@Tarny\_GG\_Channie](https://discourse.julialang.org/u/Tarny_GG_Channie)\
**Post date:** [December 1, 2023, 11:49am UTC](https://discourse.julialang.org/t/fastest-algorithm-for-xi-xj-yi-yj-for-i-in-1-n-j-in-1-n/106988/5 "2023-12-01T11:49:07Z")

</div>

Loop vectorization?

> **[GitHub - JuliaSIMD/LoopVectorization.jl: Macro(s) for vectorizing loops.](https://github.com/JuliaSIMD/LoopVectorization.jl)**
>
> Macro(s) for vectorizing loops. Contribute to JuliaSIMD/LoopVectorization.jl development by creating an account on GitHub.

---

<div class="post-metadata">

**Author:** ![heliosdrm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/heliosdrm/32/3851_2.png) [@heliosdrm](https://discourse.julialang.org/u/heliosdrm)\
**Post date:** [December 1, 2023, 12:06pm UTC](https://discourse.julialang.org/t/fastest-algorithm-for-xi-xj-yi-yj-for-i-in-1-n-j-in-1-n/106988/6 "2023-12-01T12:06:20Z")

</div>

I added `using LinearAlgebra` and changed the iterator of the inner loop by `for i in 1:j`; and then returned `Symmetric(D)` instead of just `D`. That reduced the result of `@btime` from 396.495 ms to 201.799 ms in my computer (although it did not decrease memory usage).

Two further differences: (1) your example did not return any value; (2) it should be `@btime baseline($x, $y)`, with escaped arguments, in order to report more accurate results.

---

<div class="post-metadata">

**Author:** ![aplavin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aplavin/32/222056_2.png) [@aplavin](https://discourse.julialang.org/u/aplavin)\
**Post date:** [December 1, 2023, 12:21pm UTC](https://discourse.julialang.org/t/fastest-algorithm-for-xi-xj-yi-yj-for-i-in-1-n-j-in-1-n/106988/7 "2023-12-01T12:21:02Z")

</div>

It seems already optimal, if you actually need a materialized n\*n array. There are n^2 different memory accesses and n^2 operations, so this is memory-limited. Basically, your function is already as fast as `zeros(n, n)`.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [December 1, 2023, 12:22pm UTC](https://discourse.julialang.org/t/fastest-algorithm-for-xi-xj-yi-yj-for-i-in-1-n-j-in-1-n/106988/8 "2023-12-01T12:22:06Z")

</div>

> [@juliohm](#):
>
> Is it possible to rely on FFT algorithms to compute the output efficiently?

No, this is not something an FFT can help with — you have n^2 outputs, so you can’t compute it in less than \Theta(n^2) operations.

> [@juliohm](#):
>
> Given the proposed algorithm, what is the largest `n` assuming 32Gb of RAM?

Assuming you mean 32GiB = 32 \times 2^{30} bytes, this is basically just solving 8n^2 = 32 \times 2^{30} \implies n = 2^{16} = 65536, since `sizeof(Float64) == 8`, i.e. finding the biggest matrix that fits in RAM. In practice, you can’t fill up _all_ the RAM with your matrix, since you need some memory for your OS, for Julia, etcetera. But a few tens of thousands.

Depending upon what you want to _do_ with the matrix, the best thing may be to _not store the matrix at all_, but to compute its entries as needed. For example, if you are doing linear algebra with this matrix A , then you can exploit the fact that you can compute A z (matrix–vector products) with \Theta(n) operations and \Theta(n) storage, by e.g.:

```julia
matvec(x, y, z) = (x .* y) .* sum(z) .+ sum(x .* y .* z) .- y .* (x'z) .- x .* (y'z)

```

which gives the same results as `baseline(x, y) * z` (assuming you add a `return D` to the end of your `baseline`). So, many linear-algebra calculations can be made much more efficient by using iterative algorithms based off of this.

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [December 1, 2023, 12:46pm UTC](https://discourse.julialang.org/t/fastest-algorithm-for-xi-xj-yi-yj-for-i-in-1-n-j-in-1-n/106988/9 "2023-12-01T12:46:17Z")

</div>

stevengj’s suggestion of not materializing the matrix is great. If it doesn’t work out, one imrovement might be to tile the calculation so the tiles make the calculation more cache friendly. Something like:

```julia
for j1 in 1:K:n
    for i1 in 1:K:n
        for j in j1:min(j1+K-1, n)
            for i in i1:min(i1+K-1, n)
...

```

For the right K it could make `x[i]`, `x[j]`, `y[i]`, `y[j]` accesses stay within cache.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [December 1, 2023, 12:47pm UTC](https://discourse.julialang.org/t/fastest-algorithm-for-xi-xj-yi-yj-for-i-in-1-n-j-in-1-n/106988/10 "2023-12-01T12:47:19Z")

</div>

> [@juliohm](#):
>
> that exploits parallel hardware

`@tturbo` from LoopVectorization makes a difference:

```julia-repl
julia> @btime baseline($x, $y)
  196.333 ms (2 allocations: 762.94 MiB)

julia> function baseline(x, y)
         n = length(x)
         D = Matrix{Float64}(undef, n, n)
         @tturbo for j in 1:n
           xj = x[j]
           yj = y[j]
           for i in 1:n
             xi = x[i]
             yi = y[i]
             D[i,j] = (xi - xj) * (yi - yj)
           end
         end
       end
baseline (generic function with 1 method)

julia> @btime baseline($x, $y)
  47.542 ms (2 allocations: 762.94 MiB)

julia> Threads.nthreads()
16

```

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [December 1, 2023, 2:01pm UTC](https://discourse.julialang.org/t/fastest-algorithm-for-xi-xj-yi-yj-for-i-in-1-n-j-in-1-n/106988/11 "2023-12-01T14:01:56Z")

</div>

Thank you all, what about GPUs? Can we leverage some package to convert the loops into an efficient kernel?

---

<div class="post-metadata">

**Author:** ![Tarny\_GG\_Channie](https://avatars.discourse-cdn.com/v4/letter/t/3bc359/32.png) [@Tarny\_GG\_Channie](https://discourse.julialang.org/u/Tarny_GG_Channie)\
**Post date:** [December 1, 2023, 2:06pm UTC](https://discourse.julialang.org/t/fastest-algorithm-for-xi-xj-yi-yj-for-i-in-1-n-j-in-1-n/106988/12 "2023-12-01T14:06:20Z")

</div>

Tullio?

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [December 1, 2023, 2:10pm UTC](https://discourse.julialang.org/t/fastest-algorithm-for-xi-xj-yi-yj-for-i-in-1-n-j-in-1-n/106988/13 "2023-12-01T14:10:01Z")

</div>

To give some context, this issue is showing up in a geostatistical clustering algorithm known as GHC. Given geospatial data, this algorithm constructs a dissimilarity matrix using both the geographic coordinates and the features to constraint the clusters spatially:

```julia
using GeoStats
import GLMakie

# sinusoidal field
Ω = georef((Z=[10sin(i/10) + j for i in 1:4:100, j in 1:4:100],))

Ω |> viewer

```

![image](https://global.discourse-cdn.com/julialang/original/3X/6/8/688de3307188e65bd1892cdd100cdc6a37a02a16.png)

```julia
# request 20 clusters
Ω |> GHC(20, 1.0) |> viewer

```

![image](https://global.discourse-cdn.com/julialang/original/3X/e/3/e3f332b31a33c90a413b0e6aecab1d7227799aba.png)

You can see the current implementation [here](https://github.com/JuliaEarth/GeoStatsTransforms.jl/blob/f42540b7d93cfcd2a306f6ab6bce693ca304b13e/src/clustering/ghc.jl#L105-L132), which uses `kron` products of sparse matrices to avoid blowing up memory. Any suggestion of improvement is welcome.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [December 1, 2023, 2:11pm UTC](https://discourse.julialang.org/t/fastest-algorithm-for-xi-xj-yi-yj-for-i-in-1-n-j-in-1-n/106988/14 "2023-12-01T14:11:33Z")

</div>

I will take a look at Tullio.jl, it seems great. I never had a chance to try these packages but now I have a good challenge to tackle 🙂