Outer function to identify matches: Issue with scalability

Question

I want to perform matching between two groups in a data frame consisting of 10 million rows, where all rows belonging to one group (binary) are matched with observations from the other group (with replacement) if their difference on another column is smaller than a pre-set threshold. The end result should be a data frame with 2 columns: (1) id number and (2) id number of matched row To do this, I use the outer function. See the toy example below:

set.seed(123)

# Creating data
df <- data.frame(id = c(1:10000000),
                 group = rbinom(10000000,1, 0.3),
                 value = round(runif(10000000),2))



threshold <- round(sd(df$value)*0.1,2)

#################################################################

# Identifying matches
library(tidyverse)
library(data.table)

# All values
dist_mat <- df$value

# Adding identifier
names(dist_mat) <- df$id

# Dropping combinations that are not of interest
dist_mat_col <-dist_mat[df$group == 0]
dist_mat_row <- dist_mat[df$group == 1]

# Difference between each value
dist_mat <- abs(outer(dist_mat_row, dist_mat_col, "-"))

# Identifying matches that fulfills the criteria
dist_mat <- dist_mat <= threshold 

# From matrix to a long dataframe
dist_mat <- melt(dist_mat)

# Tidying up the dataframe and dropping unneccecary columns and rows. 
dist_mat <- dist_mat %>%
  rename(id = Var1,
         matched_id = Var2,
         cond = value) %>%
  filter(cond == TRUE) %>%
  left_join(df, by = "id") %>%
  select(id, matched_id)

This code works for smaller datasets but is having issues when scaling up the data size (for obvious reasons). You can try to reduce the data frame size to 100 or 1000 rows and it should run more smoothly. The issue is related to the outer function and is stated as: Error: cannot allocate vector of size 156431.9 Gb.

As a way to solve this, I tried to do the matching row-wise, i.e., one row at a time. But this takes a tremendously long time (2500 rows in 8h, where I have 3 million rows to loop through...). See code below:

dist_mat <- df$value

names(dist_mat) <- df$id

# Dropping combinations that are not of interest
dist_mat_col <-dist_mat[df$group == 0]
dist_mat_row <- dist_mat[df$group == 1]


# Difference between each value
matched_df <- data.frame()
for (i in 1:length(dist_mat_row)) {
  print(i)
  
  
  dist_mat <- as.matrix(abs(outer(dist_mat_row[i], dist_mat_col, "-")))
  colnames(dist_mat) <- names(dist_mat_col)
  rownames(dist_mat) <- names(dist_mat_row[i])
  
  
  dist_mat <- dist_mat <= threshold 
  
  # From matrix to a long dataframe
  dist_mat <- melt(dist_mat)
  
  # Tidying up the dataframe and dropping unneccecary columns and rows. 
  dist_mat <- dist_mat %>%
    rename(id = Var1,
           matched_id = Var2,
           cond = value) %>%
    filter(cond == TRUE) %>%
    left_join(df, by = "id") %>%
    select(id, matched_id)
  
  matched_df <- rbind(matched_df, dist_mat)
  rm(dist_mat)
  gc()
}

Is there any way of doing this that does not run out of memory or takes a tremendous time? So far, I've been trying to "trim some meat" off the data to reduce the size, and perhaps there are any more ways to do this? An alternative is to not do this the "brute" way but to find an alternative. Does anyone have any suggestions or ideas?

Thanks!

Answer 1

This will be my correct answer.

First, we need a function that will generate a data set with the appropriate proportion of the number of unique values. Here it is.

library(tidyverse)
library(collapse)

fdf = function(n, nup=.1) {
  vp = 1/n/nup
  tibble(
    id = c(1:n),
    group = rbinom(n, 1, 0.3),
    value = round(runif(n)/vp)*vp)
}

For example, let's generate a set of 350 records with a ratio of unique values equal to 0.15

fdf(350, .15) %>% funique(cols=3) %>% nrow()

output

[1] 53

Now for a second example. 1000 lines with approximately 100 unique values.

fdf(1000, .1) %>% funique(cols=3) %>% nrow()

output

[1] 101

Now the most important and crucial thing. A binary search function that finds a range of val values that differ by tresh.

fbin = function(x, val, tresh = 0){
  vmin = val - tresh
  vmax = val + tresh
  n = length(x)
  e = .Machine$double.eps
  if((x[1]-vmax)>=e | (vmin-x[n])>=e) NULL else{
    l = 1
    r = n
    if(abs(x[1]-vmin)<=e | abs(x[1]-vmax)<=e |
       ((x[1]-vmin)>=e & (vmax-x[1])>=e)) imin=1 else {
      while(l <= r){
        i = (l + r)  %/% 2
        if((vmin-x[i])>e){
          l = i + 1
        } else {
          if(!(vmin-x[i-1])>e){
            r = i - 1
          } else break
        }
      }
      imin=i
    }

    l = imin
    r = n
    if(abs(x[n]-vmin)<=e | abs(x[n]-vmax)<=e |
       ((x[n]-vmin)>=e & (vmax-x[n])>=e)) imax = n else {
      while(l <= r){
        i = (l + r)  %/% 2
        if((x[i]-vmax)>e){
          r = i - 1
        } else {
          if(!((x[i+1]-vmax)>e)){
            l = l + 1
          } else break
        }
      }
      imax=i
    }
    imin:imax
  }
}

First, a few notes about this feature. I took into account the fact that the val and tresh variables of the double type, and thus, due to the inaccuracy of the calculations, ordinary comparisons cannot be used here such as x[i]>vmax or x[i]==vmax. My search function requires the argument x to be sorted in descending order!

Let's do some unit tests.

set.seed(123)
x = sample(1:10, 30, replace=T) %>% sort()
x
#[1]  1  2  3  3  3  3  3  4  4  5  5  5  6  6  7  7  7  8  9  9  9  9  9  9  9 10 10 10 10 10

x[fbin(x, 100, 0)]
#integer(0)

x[fbin(x, -10, 0)]
#integer(0)


x[fbin(x, 1, 0)]
#[1] 1

x[fbin(x, 10, 0)]
#[1] 10 10 10 10 10


x[fbin(x, 1, 1)]
#[1] 1 2

x[fbin(x, 10, 1)]
# [1]  9  9  9  9  9  9  9 10 10 10 10 10


x[fbin(x, 5, 0)]
#[1] 5 5 5

x[fbin(x, 5, 2)]
#[1] 3 3 3 3 3 4 4 5 5 5 6 6 7 7 7

x[fbin(x, 5, 10)]
# [1]  1  2  3  3  3  3  3  4  4  5  5  5  6  6  7  7  7  8  9  9  9  9  9  9  9 10 10 10 10 10

As you can see, the function returns the indexes for which the vector x values fall within the range of <val-tresh, val+tresh>.

Now it's time for a specific test. We'll see how fbin does a 10,000,000-element vector search.

set.seed(123)
n = 10000000
x = runif(n) %>% round(6) %>% sort()
funique(x) %>% length()

x[fbin(x, .5)]
#[1] 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5
x[fbin(x, .5, .000001)]
# [1] 0.499999 0.499999 0.499999 0.499999 0.499999 0.499999 0.499999 0.499999 0.499999
# [10] 0.499999 0.500000 0.500000 0.500000 0.500000 0.500000 0.500000 0.500000 0.500000
# [19] 0.500000 0.500000 0.500000 0.500000 0.500000 0.500001 0.500001 0.500001 0.500001
# [28] 0.500001 0.500001 0.500001 0.500001

Now let's see how long such a search will take.

library(microbenchmark)
ggplot2::autoplot(microbenchmark(fbin(x, .5, .001),
                                 fbin(x, .5, .002),
                                 fbin(x, .5, .003),
                                 fbin(x, .5, .004),
                                 times=10))

As you can see, the search takes about 1000 us. Now let's compare that to the subset functions.

ggplot2::autoplot(microbenchmark(x[fbin(x, .5, .001)],
                                 ss(x, x>=(0.5+0.001) & x<=(0.5-0.001)),
                                 subset(x, x>=(0.5+0.001) & x<=(0.5-0.001)),
                                 times=10))

As you can see, it is two or three orders faster!

It's time for the right function to solve your task.

fmatch = function(df, tresh){
  #Adding a column with the row number
  df = df %>% ftransform(row = 1:nrow(.))
  
  #Splitting into two sorted subsets
  df0 = df %>% roworder(value) %>% fsubset(group == 0)
  df1 = df %>% roworder(value) %>% fsubset(group == 1)
  
  #Transformations on matrices
  M0 = df0 %>% qM()
  M1 = df1 %>% qM()
  
  #Prepare unique values from group 1
  uM1 = df1$value %>% funique()
  
  out = list()
  for(i in 1:length(uM1)){
    iM0 = fbin(M0[,3], uM1[i], tresh)
    if(length(iM0)>0){
      iM1 = fbin(M1[,3], uM1[i])
      out[[paste0(uM1[i])]] = list(
        row0 = M0[iM0, 4],
        row1 = M1[iM1, 4]
      )
    }
  }
  out
}

How does this feature work? I will describe it step by step.

1. Complete the data frame with the line numbers
2. Split the frame into two sorted subsets for grooup 1 and group 0
3. Convert it into matrices to speed up the operation (maybe you don't need it)
4. Prepare unique values from the subset of groups 1
5. For each unique value in the subset of group 1, do:

5.1 In the set for group 0, search for all rows for which value does not differ from the current unique value + - threshold

5.2 If only such lines exist, write one list which will contain the line numbers from the subset for group 1 with the value equal to the current value, and the line numbers from the subset group 0.

Let's see this for an example

#Preparation of data and threshold
set.seed(123)
df = fdf(100)
threshold = round(sd(df$value)*0.1,2)

out = fmatch(df, threshold)

df[out[[1]]$row1,]
# # A tibble: 1 x 3
# id group value
# <int> <int> <dbl>
#  1    16     1   0.1

df[out[[1]]$row0,]
# # A tibble: 6 x 3
# id group value
# <int> <int> <dbl>
# 1    10     0   0.1
# 2    13     0   0.1
# 3    28     0   0.1
# 4    29     0   0.1
# 5    48     0   0.1
# 6    55     0   0.1

df[out[[2]]$row1,]
# # A tibble: 3 x 3
# id group value
# <int> <int> <dbl>
# 1    24     1   0.2
# 2    58     1   0.2
# 3    68     1   0.2

df[out[[2]]$row0,]
# # A tibble: 9 x 3
# id group value
# <int> <int> <dbl>
# 1    27     0   0.2
# 2    44     0   0.2
# 3    46     0   0.2
# 4    47     0   0.2
# 5    49     0   0.2
# 6    54     0   0.2
# 7    60     0   0.2
# 8    72     0   0.2
# 9    99     0   0.2

Now I will change the threshold to 0.2 and repeat the test.

out = fmatch(df, 0.2)

df[out[[1]]$row1,]
# # A tibble: 1 x 3
# id group value
# <int> <int> <dbl>
# 1    16     1   0.1

df[out[[1]]$row0,]
# # A tibble: 24 x 3
# id group value
# <int> <int> <dbl>
# 1    43     0   0  
# 2    10     0   0.1
# 3    13     0   0.1
# 4    28     0   0.1
# 5    29     0   0.1
# 6    48     0   0.1
# 7    55     0   0.1
# 8    27     0   0.2
# 9    44     0   0.2
# 10    46     0   0.2
# # ... with 14 more rows

df[out[[2]]$row1,]
# # A tibble: 3 x 3
# id group value
# <int> <int> <dbl>
# 1    24     1   0.2
# 2    58     1   0.2
# 3    68     1   0.2

df[out[[2]]$row0,]
# # A tibble: 32 x 3
# id group value
# <int> <int> <dbl>
# 1    43     0   0  
# 2    10     0   0.1
# 3    13     0   0.1
# 4    28     0   0.1
# 5    29     0   0.1
# 6    48     0   0.1
# 7    55     0   0.1
# 8    27     0   0.2
# 9    44     0   0.2
# 10    46     0   0.2
# # ... with 22 more rows

Now it's time to test with 100,000 rows.

set.seed(123)
df = fdf(100000)
threshold = round(sd(df$value)*0.1,2)

start_time <- Sys.time()
out = fmatch(df, threshold)
end_time <- Sys.time()

end_time - start_time
#Time difference of 13.9958 secs

object.size(out)
#319309040 bytes

As you can see, the whole thing took only 14 seconds. The output list is 320 MB. This could be crucial.

I ran another test on a set of 500,000 rows.

set.seed(123)
df = fdf(500000)
threshold = round(sd(df$value)*0.1,2)

start_time <- Sys.time()
out = fmatch(df, threshold)
end_time <- Sys.time()

end_time - start_time
#Time difference of 7.982853 mins

length(out)
#47509

object.size(out)
#7889344576 bytes

As you hang, the fivefold increase in the data set has made the time 34 times longer. The initial list has grown 24 times and now takes almost 8 GB!

There is a very important conclusion from this. Probably for 10,000,000 lines you will not have enough memory to complete the operation. So I suggest slightly modifying the fmatch function so that it returns results only for a specific subset of unique values.

Perhaps we could also optimize the binary search functionality a bit more. But I would need to know what your values are in the variable value in your dataframe.

However, as you can see, the critical factor here is not the execution time, but the memory availability.

I will be waiting for your opinion.

Also write if my solution is clear to you and if you need any additional explanations.

Last update

I did one more test tonight. However, it required minimal modification to my fmatch function. It added two additional arguments, vmin and vmax. The function will now only run for unique values in the range <vmin, vmax).

fmatch1 = function(df, tresh, vmin=0, vmax=1){
  #Adding a column with the row number
  df = df %>% ftransform(row = 1:nrow(.))

  #Splitting into two sorted subsets
  df0 = df %>% roworder(value) %>% fsubset(group == 0)
  df1 = df %>% roworder(value) %>% fsubset(group == 1)

  #Transformations on matrices
  M0 = df0 %>% qM()
  M1 = df1 %>% qM()

  #Prepare unique values from group 1
  uM1 = df1$value %>% funique() %>% ss(.>=vmin & .<vmax)

  out = list()
  for(i in 1:length(uM1)){
    iM0 = fbin(M0[,3], uM1[i], tresh)
    if(length(iM0)>0){
      iM1 = fbin(M1[,3], uM1[i])
      out[[paste0(uM1[i])]] = list(
        row0 = M0[iM0, 4],
        row1 = M1[iM1, 4]
      )
    }
  }
  out
}

Now I was able to perform a data frame test with 10,000,000 rows. However, I limited myself to values in the range <0, 0.005).

set.seed(123)
df = fdf(10000000)
threshold = round(sd(df$value)*0.1,2)

start_time <- Sys.time()
out = fmatch1(df, threshold, 0, .005)
end_time <- Sys.time()

end_time - start_time
#Time difference of 6.865934 mins

length(out)
#4706

object.size(out)
#8557061840 bytes

The whole thing took almost 7 minutes and the result was as much as 9 GB of memory !!

If we now assume that it will be relatively linear, we can expect that for all unique values in the data frame with 10,000,000 lines, the function runtime will be approx. 24 hours and the result should be approx. 1,800 GB. Unfortunately, my computer does not have that much memory.

Answer 2

In fact, what I am writing now will not be the actual answer. This is going to be quite a long comment. Unfortunately, I would not fit it in one or even several comments. Therefore, I am asking everyone to be understanding and not to criticize what I am writing here.

Now to the point. I looked at your problem. I've even been able to write a program that will do your job in much less time. With 100,000 lines, the program only ran for a few minutes. What compared to the 8 hours you gained on 2,500 rows is a clear difference. The problem, however, probably lies in the assumptions of the task itself.

When you write yourself, you have 10,000,000 rows. However, of those 10,000,000 lines, you only have 100 unique values, which is due to round(runif(n), 2)). So the first question to ask: it is the same for your real data?

Later you will say you want to match group id 0 to group id 1 if the difference between the values ​​is less than the specified threshold (let's assume the threshold for a moment is 0.3). So let's check what it gives in the output. If you only have 100 unique values ​​and 10,000,000 rows, you can expect group 0 to be around 50,000 values ​​of 0.99. Each of these values, of course, has a different id. However, in group 1, you will have approximately 3,450,000 rows with values ​​less than 0.69. Now, if you want to match each of these 50,000 IDs to 3,450,000 Group 1 IDs, you will get 172,500,000,000 matches in total !! Recall that we matched only the id from group 0, for which the value was 0.99.

Finally, my 100,000 row code generated a result set of only 10,000,000 rows! And although he did it in minutes, it strained my computer's memory a lot.

In addition, I wonder if by any chance you did not want to match the id not as you write, but when the absolute value of the difference between the values is less than the accepted threshold? abs(value1 - value0)<threshold?

If you are very curious, here is my code that I wrote about above.

library(tidyverse)
library(collapse)


set.seed(123)
n = 100000
df = tibble(
  id = c(1:n),
  group = rbinom(n,1, 0.3),
  value = round(runif(n),2))

threshold = round(sd(df$value)*0.1,2)


m1 = df %>%
  fsubset(group == 1) %>%
  roworder(value) %>%
  ftransform(row = 1:nrow(.))

m1.idx = m1 %>% funique(cols=3)
m1.M = m1 %>% qM()

m0 = df %>%
  fsubset(group == 0) %>%
  roworder(value)

m0.idx = m0 %>% funique(cols=3)
m0.M = m0 %>% qM

out = list()
for(i in 1:nrow(m0.M)){
  id0 = m0.M[i,1]
  value0 = m0.M[i,3]
  value1 = round(value0 - threshold, 2)
  idx = m1.idx %>% fsubset(value<=value1) %>% qM
  if(nrow(idx)>1){
    last.row = idx[nrow(idx), 4]-1
    out[[paste0(id0)]] = m0 %>% ss(1:last.row,1)
  }
}
dfout = unlist2d(out) %>% frename(.id = id0, id = id1) %>% qTBL()

However, I would suggest a slightly different solution. Perhaps it will be enough to remember only each of the 100 unique values from one of the groups and to each of them assign all id from group 0 for which this value exists, and all id from group 1 for which the value is less than the set threshold, or the absolute difference of values is smaller than this threshold.

Unfortunately, I do not know if such a solution would be acceptable for you. I will be waiting for a comment from you.