Speed up computation on lists of lists in R (scalability issues)
You have two lapply calls in there, so essentially two hidden for loops contributing to the slow execution.

Here, one way to tap into faster compiled functions is to rely on the data.table package for merging (a.k.a. a database join) your two data sets.

First, we turn your two lists into two data.tables:

idx1 &lt;- rep(seq_along(source.list), lengths(source.list))
idx2 &lt;- rep(seq_along(target.list), lengths(target.list))

library(data.table)
X &lt;- data.table(id = unlist(source.list), idx1 = idx1,
                pos = seq_along(idx1), key = "id")
Y &lt;- data.table(id = unlist(target.list), idx2 = idx2, key = "id")


where idx1 and idx2 are the variables telling us in which element of source.list and target.list a specific item (id) belongs. Also pos is just a temporary row number we will use later to sort the data back, as the next statement below (the merge) shuffles things around:

Z &lt;- Y[X, allow.cartesian=TRUE]
Z &lt;- Z[order(Z$pos)]


At this point, if it's possible for you, I would recommend you stop here. By that, I mean to make the next steps of your analysis use this Z data.table rather than a nested list, as to avoid further slow processing loops (for/lapply/Map/etc.). However, if you really need the nested list output, you can do:

no_name_split &lt;- function(...) unname(split(...))
gp1 &lt;- factor(Z$idx1, seq_along(source.list))
gp2 &lt;- factor(Z$idx2, seq_along(target.list))

res &lt;- Map(no_name_split, no_name_split(Z$id, gp1),
                          no_name_split(gp2,  gp1))




Below is a simulation with a larger dataset, comparing computation times and checking that the results are identical:

set.seed(632)
source.list &lt;- replicate(500, sample(100000, 1000), simplify = FALSE)
target.list &lt;- replicate(555, sample(100000, 1055), simplify = FALSE)

system.time({
  idx1 &lt;- rep(seq_along(source.list), lengths(source.list))
  idx2 &lt;- rep(seq_along(target.list), lengths(target.list))

  library(data.table)
  X &lt;- data.table(id = unlist(source.list), idx1 = idx1,
                  pos = seq_along(idx1), key = "id")
  Y &lt;- data.table(id = unlist(target.list), idx2 = idx2, key = "id")
  Z &lt;- Y[X, allow.cartesian = TRUE]
  Z &lt;- Z[order(Z$pos)]

  no_name_split &lt;- function(...) unname(split(...))
  gp1 &lt;- factor(Z$idx1, seq_along(source.list))
  gp2 &lt;- factor(Z$idx2, seq_along(target.list))
  res &lt;- Map(no_name_split, no_name_split(Z$id, gp1),
                            no_name_split(gp2,  gp1))
})
#    user  system elapsed 
#   3.394   0.382   3.646 

system.time({
  result &lt;- lapply(X = source.list, FUN = map.partition.R,
                   partitionlist = target.list)
})
#    user  system elapsed 
#  23.943   5.329  36.999 

identical(res, result)
# [1] TRUE