# Efficient way of computing norm between all the rows of a matrix

**URL:** https://discourse.julialang.org/t/efficient-way-of-computing-norm-between-all-the-rows-of-a-matrix/70277
**Category:** New to Julia
**Tags:** question, linearalgebra
**Created:** [October 24, 2021, 12:10am UTC](https://discourse.julialang.org/t/efficient-way-of-computing-norm-between-all-the-rows-of-a-matrix/70277 "2021-10-24T00:10:39Z")
**Posts on this page:** 6
**Page:** 1

<div class="post-metadata">

### Author: ![vinicius\_de\_lima](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vinicius_de_lima/32/28717_2.png) [@vinicius\_de\_lima](https://discourse.julialang.org/u/vinicius_de_lima)
#### Post date: [October 24, 2021, 12:10am UTC](https://discourse.julialang.org/t/efficient-way-of-computing-norm-between-all-the-rows-of-a-matrix/70277/1 "2021-10-24T00:10:39Z")

</div>

I stored points in the rows of a `999900 x 2` matrix and tried to count the number of pairs of these points that are at a distance smaller than `r`, but my implementation was too slow, there is a more efficient way to achieve this? Below is my attempt:

```julia
function C(pontos::T,r::Real) where T <: AbstractMatrix{<:Real}
    contar = 0
    for r1 in 1:size(pontos,1)
        for r2 in r1+1:size(pontos,1)
            if norm(pontos[r1,:] - pontos[r2,:]) < r
                contar += 1 
            end
        end
    end
    return contar
end

```

---

<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: [October 24, 2021, 1:06am UTC](https://discourse.julialang.org/t/efficient-way-of-computing-norm-between-all-the-rows-of-a-matrix/70277/2 "2021-10-24T01:06:08Z")

</div>

Try changing your data structure to a 1d array of 2-component SVectors (from StaticArrays.jl). That is a vastly more efficient way to do computations with small fixed-length vectors (e.g. coordinates in 2d).

For example, it will eliminate all the heap allocations you are doing in your inner loops, and it will inline all of the 2-component vector operations (unrolling the 2-iteration loops).

---

<div class="post-metadata">

### Author: ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [October 24, 2021, 2:06am UTC](https://discourse.julialang.org/t/efficient-way-of-computing-norm-between-all-the-rows-of-a-matrix/70277/3 "2021-10-24T02:06:23Z")

</div>

also, computing norm squared is generally faster since it doesn’t require a square root.

---

<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: [October 24, 2021, 2:50am UTC](https://discourse.julialang.org/t/efficient-way-of-computing-norm-between-all-the-rows-of-a-matrix/70277/4 "2021-10-24T02:50:20Z")

</div>

If `r` is significantly smaller than the average distance of the points, you can do that much, much faster using cell lists. For instance using [CellListMap.jl](https://github.com/m3g/CellListMap.jl). if you can provide an actual example and need some help, let me know.

---

<div class="post-metadata">

### Author: ![vinicius\_de\_lima](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vinicius_de_lima/32/28717_2.png) [@vinicius\_de\_lima](https://discourse.julialang.org/u/vinicius_de_lima)
#### Post date: [October 24, 2021, 4:49am UTC](https://discourse.julialang.org/t/efficient-way-of-computing-norm-between-all-the-rows-of-a-matrix/70277/5 "2021-10-24T04:49:01Z")

</div>

Thank you, I believe that in fact `r` will attend this condition, in my use case it will assume values in the range `2^-8` to `2^-2` and the points are from a [Henon map](https://en.wikipedia.org/wiki/H%C3%A9non_map), [here](https://github.com/vinicius-souza-lima/FM-223-Dinamica-Nao-Linear-e-Caos/blob/master/Lista%2006/pontos%20mapa%20de%20henon.txt) is the matrix I used.

---

<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: [October 24, 2021, 12:28pm UTC](https://discourse.julialang.org/t/efficient-way-of-computing-norm-between-all-the-rows-of-a-matrix/70277/6 "2021-10-24T12:28:44Z")

</div>

This kind of data is not the one where cell lists shine the most, because it is very sparse (cell lists are ideal for more or less homogeneous distributions of points). Still, there it goes (I have adapted the naive version to use static arrays):

```julia
julia> include("./vinicius.jl")
cell list - serial: 
  5.258 s (13576 allocations: 68.23 MiB)
cell list - parallel: 
  2.096 s (77501 allocations: 131.05 MiB)
naive - serial: 
# still running, will update when finished... 

```

Here is the code:

```julia
using CellListMap
using StaticArrays
using DelimitedFiles
using LinearAlgebra: norm
using BenchmarkTools

function count_close(p::AbstractVector{SVector{N,T}},r;parallel=true) where {N,T}
    p_length = limits(p) # length of set on each dimension
    # To prevent too many cells, set the cell side such that the maximum number of cells 
    # is at most equal to the number of points
    max_length = maximum(p_length.limits) 
    max_bins = max_length/r
    n = ceil(Int,min(max_bins,length(p)^(1/N))) 
    cell_side = max(r,max_length/n)
    # Prepare cell lists
    box = Box(p_length,cell_side)
    cl = CellList(p,box,parallel=parallel)
    r2 = r^2
    n_close = 0
    n_close = map_pairwise!(
        (x,y,i,j,d2,n_close) -> d2 < r2 ? n_close + 1 : n_close,
        n_close,
        box, cl, parallel=parallel
    )
    return n_close
end

data = readdlm("pontos.txt")
pontos = [SVector{2,Float64}(data[i,1],data[i,2]) for i in axes(data,1) ]
pontos = pontos[1:100000]
r = 2e-2

println("cell list - serial: ")
@btime count_close($pontos,$r,parallel=false)

println("cell list - parallel: ")
@btime count_close($pontos,$r,parallel=true)

function C(pontos,r)
    contar = 0
    for r1 in 1:length(pontos)
        for r2 in r1+1:length(pontos)
            if norm(pontos[r2]-pontos[r1]) < r
                contar += 1 
            end
        end
    end
    return contar
end

println("naive - serial: ")
@btime C(pontos,r)

```

up: For 100\_000 points this is the difference (for one million I think it will never end, had to kill it):

```julia
julia> include("./vinicius.jl")
cell list - serial: 
  74.476 ms (8778 allocations: 7.61 MiB)
cell list - parallel: 
  27.496 ms (40732 allocations: 17.19 MiB)
naive - serial: 
  8.912 s (1 allocation: 16 bytes)
15579644

```
