Pairwise Differences

This notebook is based on R code in a blog post by Erwin Kalvelagen: “Arranging points on a line”. The problem considered in that post involves placing points on a line so as to maximize the minimum difference between any pair of consecutive points. Our focus here is not on the problem itself, but on the objective function for a GA coded in R.

We will compare two computational approaches, looking for differences in execution time. The first approach, using Erwin’s code, is to compute absolute values of all pairwise differences using nested loops. The second is to sort the candidate solution and then compute differences only for consecutive values.

To time the methods, we will use the tictoc package.

library(tictoc)

The argument to each of the following functions consists of an unsorted vector of values on the real number line. The nested loop approach uses code from Erwin’s post.

obj <- function(x) {
  smallest <- Inf
  n <- length(x)
  for(i in 1:n)
    if (i<n)
      for(j in (i+1):n) {
        smallest <- min(abs(x[i]-x[j]),smallest)
      }
  smallest
}

The alternative approach is as follows. (Note that the |> pipe operator is new in R 4.1.0. It can be replaced by the pipe operator %>% from the magrittr package, ), or by nesting the various function calls.

obj2 <- function(x) {
  x |> sort() |> diff() |> min()
}

As a sanity check, we compare the result using a random vector of dimension 50.

x <- runif(n = 50, min = -10, max = 10)
obj(x)

[1] 0.02013424

obj2(x)

[1] 0.02013424

The results agree.

We now define a function that generates a random vector of given dimension, computes the objective value using both functions, and prints their respective execution times. We assign the result to a dummy variable to avoid having the time taken to print it count toward the method.

test <- function(n) {
  x <- runif(n, min = -10, max = 10)
  tic("Looping")
  y1 <- obj(x)
  toc()
  tic("Sorting")
  y2 <- obj2(x)
  toc()
  # Make sure the results agree.
  if (y1 != y2) {
    cat("Oops: different results ", y1, " and ", y2, "!\n")
  }
}

We try vectors of various dimensions.

for (n in c(50, 500, 5000, 10000)) {
  cat("\nVector dimension = ", n, ":\n")
  test(n)
}

Vector dimension =  50 :
Looping: 0.001 sec elapsed
Sorting: 0.003 sec elapsed

Vector dimension =  500 :
Looping: 0.065 sec elapsed
Sorting: 0 sec elapsed

Vector dimension =  5000 :
Looping: 5.966 sec elapsed
Sorting: 0 sec elapsed

Vector dimension =  10000 :
Looping: 23.268 sec elapsed
Sorting: 0.001 sec elapsed

The time required for the looping approach is quadratic in the vector dimension. The time for the sorting approach is so short it’s hard to believe. We can try higher dimensions with just the sorting approach to get a bit better idea of its speed.

test2 <- function(n) {
  x <- runif(n, min = -10, max = 10)
  tic("Sorting")
  y2 <- obj2(x)
  toc()
}

for (n in c(50000, 100000, 500000, 1000000, 2000000)) {
  test2(n)
}

Sorting: 0.003 sec elapsed
Sorting: 0.007 sec elapsed
Sorting: 0.046 sec elapsed
Sorting: 0.113 sec elapsed
Sorting: 0.332 sec elapsed

With a vector of 2,000,000 observations, evaluation time is still well under one second.

There is one remaining question. The biggest difference in execution time stems from the looping approach being \(O(n^2)\) and the sorting approach being \(O(n \log n)\). That said, some people contend that explicit loops in R are significantly less efficient that using vectorized functions (which will still loop, but with loops presumably coded in C or C++ in the function definitions). We can test that to some extent using yet another version of the objective function, which does not sort and computes all pairwise differences, but using vectorized methods.

obj3 <- function(x) {
  # Compute all pairwise absolute differences.
  d <- outer(x, x, FUN = function(y, z) abs(y - z))
  # Eliminate the zero diagonal entries.
  diag(d) <- Inf
  # Return the minimum value.
  min(d)
}

Before proceeding, we confirm on our test vector that the result is correct.

obj(x)

[1] 0.02013424

obj2(x)

[1] 0.02013424

obj3(x)

[1] 0.02013424

We can now time all three methods on random test vectors.

test3 <- function(n) {
  x <- runif(n, min = -10, max = 10)
  tic("Looping")
  y1 <- obj(x)
  toc()
  tic("Sorting")
  y2 <- obj2(x)
  toc()
  tic("Outer product")
  y3 <- obj3(x)
  toc()
  # Make sure the results agree.
  if (y1 != y3 || y2 != y3) {
    cat("Oops: different results ", y1, ", ", y2, " and ", y3, "!\n")
  }
}

for (n in c(50, 500, 5000, 10000)) {
  cat("\nVector dimension = ", n, ":\n")
  test3(n)
}

Vector dimension =  50 :
Looping: 0.002 sec elapsed
Sorting: 0.001 sec elapsed
Outer product: 0 sec elapsed

Vector dimension =  500 :
Looping: 0.062 sec elapsed
Sorting: 0.001 sec elapsed
Outer product: 0.002 sec elapsed

Vector dimension =  5000 :
Looping: 5.525 sec elapsed
Sorting: 0 sec elapsed
Outer product: 0.667 sec elapsed

Vector dimension =  10000 :
Looping: 21.413 sec elapsed
Sorting: 0.001 sec elapsed
Outer product: 2.149 sec elapsed

The outer product approach is about a factor of 10 faster than the nested loops (perhaps less if we go to higher dimensions), even though the outer product calculates every difference twice (and the difference of each observation with itself), doing slightly more than double the work of the nested loop. So there is some merit to the claim that for loops in R are inefficient.