Alternative to computation in a loop for distances for fast performance

Some ways to do this are:

pairwise2(x::AbstractMatrix, y::AbstractMatrix) = sum(abs2, x; dims=1)' .+ sum(abs2, y; dims=1) .- 2 .* (x'*y)  # written to use columns
Dist ≈ pairwise2(Points[order, :]', newP')

using Distances
Dist ≈ pairwise(SqEuclidean(), Points[order, :], newP; dims=1)  # dims=1 means rows here

using Tullio
Dist ≈ @tullio dist_t[i] := (Points[order[i], k] - newP[k])^2

Times:

julia> using BenchmarkTools

julia> @btime for i in eachindex($order)
               $Dist[i] = sum(($Points[$order[i], :]' - $newP) .^ 2)
       end
  min 514.084 μs, mean 577.836 μs (20100 allocations, 1.53 MiB)

julia> @btime pairwise2($Points[$order, :]', $newP');
  min 31.833 μs, mean 64.153 μs (11 allocations, 471.55 KiB)

julia> @btime pairwise(SqEuclidean(), $Points[$order, :], $newP; dims=1)
  min 27.459 μs, mean 45.095 μs (7 allocations, 262.02 KiB)

julia> @btime @tullio dist_t[i] := ($Points[$order[i], k] - $newP[k])^2;
  min 20.000 μs, mean 24.760 μs (2 allocations, 52.42 KiB)

julia> using LoopVectorization

julia> @btime @tullio dist_t[i] := ($Points[$order[i], k] - $newP[k])^2;
  min 13.375 μs, mean 17.343 μs (2 allocations, 52.42 KiB)
7 Likes