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)