find neighbouring elements of a matrix in R

Question

Edit: huge thanks to the users below for great contributions and to Gregor for benchmarking.

Say I have a matrix filled with integer values like this...

    mat <- matrix(1:100, 10, 10)

I can create a list of x, y coordinates of each element like this...

    addresses <- expand.grid(x = 1:10, y = 1:10)

Now for each of these coordinates (i.e. for each element in mat) i would like to find the neighboring elements (including diagonals this should make 8 neighbors).

I'm sure there is an easy way, can anyone help?

What I have tried so far is to loop through and for each element record the neighboring elements as follows;

    neighbours <- list()
    for(i in 1:dim(addresses)[1]){
      x <- addresses$x[i]
      y <- addresses$y[i]
      neighbours[[i]] <- c(mat[y-1, x  ],
                           mat[y-1, x+1],
                           mat[y  , x+1],
                           mat[y+1, x+1],
                           mat[y+1, x  ],
                           mat[y+1, x-1],
                           mat[y  , x-1],
                           mat[y-1, x-1])
    }

This runs into problems when it hits the edge of the matrix, particularly when the index is greater than the edge of the matrix.

What do you want to do at the edges? You could "wrap" and consider the left edge adjacent to the right edge, or your could return `NA`, or maybe something else. — Gregor Thomas, Mar 17 '15 at 17:14
If the answer is `NA`, you could avoid a lot of special handling in the code if you pad your matrix with `NA`s (make it 12x12 with the first and last row and column all NA), then you could just look at all neighbors of the "middle". — Gregor Thomas, Mar 17 '15 at 17:16

Gregor Thomas · Answer 1 · 2015-03-17T20:57:48.967

Here's a nice example. I did 4x4 so we can see it easily, but it's all adjustable by n. It's also fully vectorized so should have good speed.

n = 4
mat = matrix(1:n^2, nrow = n)
mat.pad = rbind(NA, cbind(NA, mat, NA), NA)

With the padded matrix, the neighbors are just n by n submatrices, shifting around. Using compass directions as labels:

ind = 2:(n + 1) # row/column indices of the "middle"
neigh = rbind(N  = as.vector(mat.pad[ind - 1, ind    ]),
              NE = as.vector(mat.pad[ind - 1, ind + 1]),
              E  = as.vector(mat.pad[ind    , ind + 1]),
              SE = as.vector(mat.pad[ind + 1, ind + 1]),
              S  = as.vector(mat.pad[ind + 1, ind    ]),
              SW = as.vector(mat.pad[ind + 1, ind - 1]),
              W  = as.vector(mat.pad[ind    , ind - 1]),
              NW = as.vector(mat.pad[ind - 1, ind - 1]))

mat
#      [,1] [,2] [,3] [,4]
# [1,]    1    5    9   13
# [2,]    2    6   10   14
# [3,]    3    7   11   15
# [4,]    4    8   12   16

  neigh[, 1:6]
#    [,1] [,2] [,3] [,4] [,5] [,6]
# N    NA    1    2    3   NA    5
# NE   NA    5    6    7   NA    9
# E     5    6    7    8    9   10
# SE    6    7    8   NA   10   11
# S     2    3    4   NA    6    7
# SW   NA   NA   NA   NA    2    3
# W    NA   NA   NA   NA    1    2
# NW   NA   NA   NA   NA   NA    1

So you can see for the first element mat[1,1], starting at North and going clockwise the neighbors are the first column of neigh. The next element is mat[2,1], and so on down the columns of mat. (You can also compare to @mrip's answer and see that our columns have the same elements, just in a different order.)

Benchmarking

Small matrix

mat = matrix(1:16, nrow = 4)
mbm(gregor(mat), mrip(mat), marat(mat), u20650(mat), times = 100)
# Unit: microseconds
#         expr     min       lq      mean   median       uq      max neval  cld
#  gregor(mat)  25.054  30.0345  34.04585  31.9960  34.7130   61.879   100 a   
#    mrip(mat) 420.167 443.7120 482.44136 466.1995 483.4045 1820.121   100   c 
#   marat(mat) 746.462 784.0410 812.10347 808.1880 832.4870  911.570   100    d
#  u20650(mat) 186.843 206.4620 220.07242 217.3285 230.7605  269.850   100  b

On a larger matrix I had to take out user20650's function because it tried to allocate a 232.8 Gb vector, and I also took out Marat's answer after waiting for about 10 minutes.

mat = matrix(1:500^2, nrow = 500)

mbm(gregor(mat), mrip(mat), times = 100)
# Unit: milliseconds
#         expr       min        lq      mean    median        uq      max neval cld
#  gregor(mat) 19.583951 21.127883 30.674130 21.656866 22.433661 127.2279   100   b
#    mrip(mat)  2.213725  2.368421  8.957648  2.758102  2.958677 104.9983   100  a

So it looks like in any case where time matters, @mrip's solutions is by far the fastest.

Functions used:

gregor = function(mat) {
    n = nrow(mat)
    mat.pad = rbind(NA, cbind(NA, mat, NA), NA)
    ind = 2:(n + 1) # row/column indices of the "middle"
    neigh = rbind(N  = as.vector(mat.pad[ind - 1, ind    ]),
                  NE = as.vector(mat.pad[ind - 1, ind + 1]),
                  E  = as.vector(mat.pad[ind    , ind + 1]),
                  SE = as.vector(mat.pad[ind + 1, ind + 1]),
                  S  = as.vector(mat.pad[ind + 1, ind    ]),
                  SW = as.vector(mat.pad[ind + 1, ind - 1]),
                  W  = as.vector(mat.pad[ind    , ind - 1]),
                  NW = as.vector(mat.pad[ind - 1, ind - 1]))
    return(neigh)
}

mrip = function(mat) {
    m2<-cbind(NA,rbind(NA,mat,NA),NA)
    addresses <- expand.grid(x = 1:4, y = 1:4)
    ret <- c()
    for(i in 1:-1)
        for(j in 1:-1)
            if(i!=0 || j !=0)
                ret <- rbind(ret,m2[addresses$x+i+1+nrow(m2)*(addresses$y+j)]) 
    return(ret)
}

get.neighbors <- function(rw, z, mat) {
    # Convert to absolute addresses 
    z2 <- t(z + unlist(rw))
    # Choose those with indices within mat 
    b.good <- rowSums(z2 > 0)==2  &  z2[,1] <= nrow(mat)  &  z2[,2] <= ncol(mat)
    mat[z2[b.good,]]
}

marat = function(mat) {
    n.row = n.col = nrow(mat)
    addresses <- expand.grid(x = 1:n.row, y = 1:n.col)
    # Relative addresses
    z <- rbind(c(-1,0,1,-1,1,-1,0,1), c(-1,-1,-1,0,0,1,1,1))
    apply(addresses, 1,
          get.neighbors, z = z, mat = mat) # Returns a list with neighbors
}

u20650 = function(mat) {
    w <-  which(mat==mat, arr.ind=TRUE)
    d <- as.matrix(dist(w, "maximum", diag=TRUE, upper=TRUE))
    # extract neighbouring values for each element
    # extract where max distance is one
    a <- apply(d, 1, function(i) mat[i == 1] )
    names(a)  <- mat
    return(a)
}

Interesting, and good analysis. I wonder why mine is faster for large matrices, I would have thought yours would be because there is only a single call to `rbind`, which should mean less allocating and copying. — mrip, Mar 18 '15 at 16:41
I was surprised too. Maybe the `as.vector` conversion is slow? — Gregor Thomas, Mar 18 '15 at 16:53
Great answer and great job, but I think you forgot to change `expand.grid` in the `mrip` function; for this reason it is the fastest. `gregor` function is actually much faster. — Leonardo, Jan 22 '21 at 14:43

score 9 · Answer 2 · answered Mar 17 '15 at 17:34

This will get you a matrix with columns corresponding to neighbors of each entry in the matrix:

mat <- matrix(1:16, 4, 4)
m2<-cbind(NA,rbind(NA,mat,NA),NA)
addresses <- expand.grid(x = 1:4, y = 1:4)

ret<-c()
for(i in 1:-1)
  for(j in 1:-1)
    if(i!=0 || j !=0)
      ret<-rbind(ret,m2[addresses$x+i+1+nrow(m2)*(addresses$y+j)]) 


> ret
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
[1,]    6    7    8   NA   10   11   12   NA   14    15    16    NA    NA    NA
[2,]    2    3    4   NA    6    7    8   NA   10    11    12    NA    14    15
[3,]   NA   NA   NA   NA    2    3    4   NA    6     7     8    NA    10    11
[4,]    5    6    7    8    9   10   11   12   13    14    15    16    NA    NA
[5,]   NA   NA   NA   NA    1    2    3    4    5     6     7     8     9    10
[6,]   NA    5    6    7   NA    9   10   11   NA    13    14    15    NA    NA
[7,]   NA    1    2    3   NA    5    6    7   NA     9    10    11    NA    13
[8,]   NA   NA   NA   NA   NA    1    2    3   NA     5     6     7    NA     9
     [,15] [,16]
[1,]    NA    NA
[2,]    16    NA
[3,]    12    NA
[4,]    NA    NA
[5,]    11    12
[6,]    NA    NA
[7,]    14    15
[8,]    10    11

how could you tweak this to get n-depth neighbors. for example to increase `ret` to also include second, third etc nearest neighbors? — Rafael, Dec 02 '19 at 17:58

user20650 · Answer 3 · 2015-03-17T20:08:14.137

Perhaps you may be able to use a distance function here using the row and column indices of the matrix elements.

# data
(mat <- matrix(16:31, 4, 4))
     [,1] [,2] [,3] [,4]
[1,]   16   20   24   28
[2,]   17   21   25   29
[3,]   18   22   26   30
[4,]   19   23   27   31

# find distances between row and column indexes
# interested in values where the distance is one
w <-  which(mat==mat, arr.ind=TRUE)
d <- as.matrix(dist(w, "maximum", diag=TRUE, upper=TRUE))

# extract neighbouring values for each element
# extract where max distance is one
a <- apply(d, 1, function(i) mat[i == 1] )
names(a)  <- mat
a

$`16`
[1] "17" "20" "21"

$`17`
[1] "16" "18" "20" "21" "22"

$`18`
[1] "17" "19" "21" "22" "23
... ....
... ....

Needs tidied, but maybe gives an idea

score 1 · Answer 4 · answered Mar 17 '15 at 18:56

Here is another approach:

n.col <- 5
n.row <- 10
mat <- matrix(seq(n.col * n.row), n.row, n.col)

addresses <- expand.grid(x = 1:n.row, y = 1:n.col)

# Relative addresses
z <- rbind(c(-1,0,1,-1,1,-1,0,1),c(-1,-1,-1,0,0,1,1,1))

get.neighbors <- function(rw) {
  # Convert to absolute addresses 
  z2 <- t(z + unlist(rw))
  # Choose those with indices within mat 
  b.good <- rowSums(z2 > 0)==2  &  z2[,1] <= nrow(mat)  &  z2[,2] <=ncol(mat)
  mat[z2[b.good,]]
}

apply(addresses,1, get.neighbors) # Returns a list with neighbors

score 1 · Answer 5 · answered Dec 27 '21 at 22:24

We can try the following base R code

f <- Vectorize(function(mat, x, y) {
  X <- pmin(pmax(x + c(-1:1), 1), nrow(mat))
  Y <- pmin(pmax(y + c(-1:1), 1), ncol(mat))
  mat[as.matrix(subset(
    unique(expand.grid(X = X, Y = Y)),
    !(X == x & Y == y)
  ))]
},vectorize.args = c("x","y"))

neighbors <- with(addresses, f(mat,x,y))

and we will see

> head(neighbors)
[[1]]
[1]  2 11 12

[[2]]
[1]  1  3 11 12 13

[[3]]
[1]  2  4 12 13 14

[[4]]
[1]  3  5 13 14 15

[[5]]
[1]  4  6 14 15 16

[[6]]
[1]  5  7 15 16 17

score 0 · Answer 6 · answered Dec 27 '21 at 22:01

There was also an much faster solution posted here: List of n first Neighbors from a 3d Array R

for large matrices this solution is much faster, see my benchmark here:

For small matrices I want to show another very interesting solution, which is fast for small matrices, but slow for bigger ones:

Using complex numbers and the abs function.

get_neighbor <- function(matrix, x=1,y=1){

    z <- complex(real = rep(1:nrow(matrix), ncol(matrix)),
                 imaginary = rep(1:ncol(matrix), each = nrow(matrix)))

    lookup <- lapply(seq_along(z), function(x){
      ## calculate distance
      dist <- which(abs(z - z[x]) < 2)
      ## remove those with dist == 0 -> it´s the number itself
      dist[which(dist != x)]
    })
    index <- (y-1)*(nrow(matrix))+x
    matrix[lookup[[index]]]

}

find neighbouring elements of a matrix in R

6 Answers6

Benchmarking

Linked