How to vectorize a function

Viewed 178

I have a 5x4 matrix. I have created a function call fun1, fun1 use double for loop to loop through the matrix and use distance function to work out the distance between two-row. The final results matrix will be a 5x5 matrix. I am struggling to covert this fun1 to a vectorization function(no loop, only apply function).

x = 
     [,1] [,2] [,3] [,4]
[1,]    1    6   11   16
[2,]    2    7   12   17
[3,]    3    8   13   18
[4,]    4    9   14   19
[5,]    5   10   15   20
distance = function(a, b) {
  sqrt(sum((a - b)^2))
}

fun1 = function(x) {
  n = nrow(x)
  results = matrix(0, nrow = n, ncol = n)
  for (i in seq_len(n)) {
    for (j in seq_len(n)) {
      results[i,j] = distance(m[i,], m[j,])
    }    
  }
  results
}
2 Answers

You can do it with just a matrix multiplication, some additions and a transpose.

x <- matrix(1:20, nrow = 5)

z <- x %*% t(x)
sqrt(diag(z)+t(diag(z)-2*z))

#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]    0    2    4    6    8
#> [2,]    2    0    2    4    6
#> [3,]    4    2    0    2    4
#> [4,]    6    4    2    0    2
#> [5,]    8    6    4    2    0

Interestingly this is faster than the in built method mentioned in the comments above!

mdist <- function(x) {
  z <- x %*% t(x)
  sqrt(diag(z)+t(diag(z)-2*z))
}

n <- 1000
l <- 100
x <- matrix(runif(n*l), ncol = l)
microbenchmark::microbenchmark(
  z1 = as.matrix(dist(x)),
  z2 = dist(x, diag = TRUE, upper = TRUE),
  z3 = mdist(x),
  times = 100
)
#> Unit: milliseconds
#>  expr      min       lq     mean   median        uq      max neval
#>    z1 82.98502 90.20049 98.54552 94.85027 101.78114 140.1809   100
#>    z2 72.54279 76.22054 82.75410 79.31865  83.47765 231.3008   100
#>    z3 54.58258 59.73461 65.62313 63.14435  67.49865 115.0379   100

In a pinch, Vectorize can do what you need:

outer(seq_len(nrow(m)), seq_len(nrow(m)),
      Vectorize(function(i,j) distance(m[i,], m[j,]), vectorize.args=c("i","j")))
#      [,1] [,2] [,3] [,4] [,5]
# [1,]    0    2    4    6    8
# [2,]    2    0    2    4    6
# [3,]    4    2    0    2    4
# [4,]    6    4    2    0    2
# [5,]    8    6    4    2    0

Vectorize takes a function as an argument and returns a function that accepts vectors, iterating internally. The function passed to it is called once for each element within the vector passed. By default, Vectorize only vectorizes the first argument of the function, but it can "zip" along multiple arguments, assuming they are all the same length, by using vectorize.args=.

This might be a little easier to visualize by redefining distance:

distance_ind = function(i, j, data) {
  sqrt(sum((data[i,] - data[j,])^2))
}
distance_ind(1, 2, m)
# [1] 2

distance_ind(c(1,3), c(2,3), m)
# [1] 2 ### wrong
distance_ind_vec <- Vectorize(distance_ind, vectorize.args = c("i", "j"))
distance_ind_vec(c(1,3), c(2,3), m)
# [1] 2 0

And the outer call:

outer(seq_len(nrow(m)), seq_len(nrow(m)), distance_ind_vec, data = m)
#      [,1] [,2] [,3] [,4] [,5]
# [1,]    0    2    4    6    8
# [2,]    2    0    2    4    6
# [3,]    4    2    0    2    4
# [4,]    6    4    2    0    2
# [5,]    8    6    4    2    0
Related