Adjust the max load factor of the dictionary to 50%

[DavisVaughan]

I've been studying our dictionary implementation to see if I can improve it in any way, since it is used in so many functions.

One thing that stood out was that our max load factor is 77%. I think this is quite high, actually. Base R uses 50%, FWIW https://github.com/wch/r-source/blob/e3f31ad3039eaade6bf9bb0d9201eb13ea33cf66/src/main/unique.c#L439.

I found a few images that show that as the load factor gets close to ~70%, the number of collisions sharply increases. The medium article suggests that quadratic probing (which we use) has a surprisingly bad result for "average time to figure out if the table contains x" as the load factor increases above 50%.
https://cs.stackexchange.com/questions/10273/hash-table-collision-probability
https://medium.com/geekculture/performance-of-hash-implementations-on-various-workloads-fedac579a39b#f31c

With all this in mind, I did some experiments with a 50% load factor and was pretty happy with what I saw. I suggest that we move to either a 50% or 60% load factor unconditionally. I think anything less than ~65% is probably going to result in significant improvement, so I'm ok if we don't want to move all the way to 50%, but I don't see much downside in doing so.

I'll show a number of examples below. I temporarily added a little helper that lets me dynamically switch the load factor with with_load_50() and with_load_77() to easily show the differences. You can see it in the commits.

---

This is particularly meaningful for functions that look up values in one vector using a dictionary built on another vector (like vec_match/in() and vec_set_intersect()). The dictionary is built on haystack for vec_match(), and it is completely reasonable for each element of haystack to be unique, so as the size of the haystack approaches the 77% load factor size, that means there is only 23% of the table left for "new" values in needles to get mapped to, which results in a lot of collisions as values in needles are first mapped to an existing haystack value, compared against it to find it is not the same value, and then quadratic probing generates a new probe, and this continues until we hit one of those empty slots in the 23% of the table that is still free.

In this first example, I'll use a haystack of size 1,600,000 with all unique values. With our current approach, this generates:

# Dict key size
next_power_of_2(1600000 / .77) = 2097152
# Load factor
1600000 / 2097152 = 76.3%

If you do the same computations with a 50% max load factor, you get:

# Dict key size
next_power_of_2(1600000 / .5) = 4194304
# Load factor
1600000 / 4194304 = 38.1%

The needles sample from 10,000,000 values, so a much wider range than haystack, meaning there will be a lot of unmatched values.

# integers
needles <- sample(10000000, 1000000, T)
haystack <- sample(1600000, 1600000, F)

bench::mark(
  "0.50" = with_load_50(vec_match(needles, haystack)),
  "0.77" = with_load_77(vec_match(needles, haystack)),
  iterations = 100
)
#> # A tibble: 2 × 6
#>   expression      min   median `itr/sec` mem_alloc `gc/sec`
#>   <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
#> 1 0.50         84.9ms   90.7ms     10.9     29.7MB    11.4 
#> 2 0.77        132.4ms  155.8ms      6.49    21.7MB     6.23

# doubles
needles <- needles + 0
haystack <- haystack + 0

bench::mark(
  "0.50" = with_load_50(vec_match(needles, haystack)),
  "0.77" = with_load_77(vec_match(needles, haystack)),
  iterations = 100
)
#> # A tibble: 2 × 6
#>   expression      min   median `itr/sec` mem_alloc `gc/sec`
#>   <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
#> 1 0.50          113ms    120ms      8.40    29.7MB    61.6 
#> 2 0.77          234ms    246ms      4.06    21.7MB     3.46

This also extends out to data frames. For data frames the equality check is more expensive, so if we can reduce the number of collisions that directly reduces the number of equality checks, which is a double win! This data frame example again shows the "worst" case of all unique values + close to 77% load factor, and then it also shows a slightly less bad case where the haystack can have duplicates, lowering the current approach's load factor to around 63%, but the new approach still does better.

set.seed(123)

# 4 mil possible combinations
dict <- tidyr::expand_grid(
  x = seq_len(2000),
  y = seq_len(2000)
)

# needles = 10 mil, with replacement
# haystack = 1.6 mil, all unique, so same 76% load factor
needles <- dplyr::sample_n(dict, size = 10000000, replace = TRUE)
haystack <- dplyr::sample_n(dict, size = 1600000, replace = FALSE)

bench::mark(
  "0.50" = with_load_50(vec_match(needles, haystack)),
  "0.77" = with_load_77(vec_match(needles, haystack)),
  iterations = 10
)
#> # A tibble: 2 × 6
#>   expression      min   median `itr/sec` mem_alloc `gc/sec`
#>   <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
#> 1 0.50        956.8ms  977.7ms     1.02     98.4MB   0.255 
#> 2 0.77           1.4s    1.51s     0.657    90.4MB   0.0730

# haystack = 1.6 mil, with replacement, so lower load factor
haystack <- dplyr::sample_n(dict, size = 1600000, replace = TRUE)

# 2 ^ 20 = 1048576
# 2 ^ 21 = 2097152
# 2 ^ 22 = 4194304
vec_unique_count(haystack) / .77
#> [1] 1711457
vec_unique_count(haystack) / .50
#> [1] 2635644

# 63% (max 77%) vs 31% (max 50%) load factors
vec_unique_count(haystack) / 2097152
#> [1] 0.6283865
vec_unique_count(haystack) / 4194304
#> [1] 0.3141932

bench::mark(
  "0.50" = with_load_50(vec_match(needles, haystack)),
  "0.77" = with_load_77(vec_match(needles, haystack)),
  iterations = 10
)
#> # A tibble: 2 × 6
#>   expression      min   median `itr/sec` mem_alloc `gc/sec`
#>   <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
#> 1 0.50       845.77ms 870.35ms     1.15     98.4MB    0.289
#> 2 0.77          1.17s    1.23s     0.819    90.4MB    0.205

---

Now lets take a look at an example where in theory our changes shouldn't help much at all.

With 9,500,000 unique haystack values, we'd get:

# 77%:
# Dict key size
next_power_of_2(9500000 / .77) = 16777216
# Load factor
9500000 / 16777216 = 56.6%

# 50%:
# Dict key size
next_power_of_2(9500000 / .50) = 33554432
# Load factor
9500000 / 33554432 = 28.3%

So the load factor would hit 56% in the current approach, which really shouldn't hurt performance much. And we indeed do see that the timings aren't very different here. The 50% load does slightly better, at the cost of more memory because the next power of 2 is much larger so we need a larger dictionary.

set.seed(1234)

needles <- sample(10000000, 1000000, replace = TRUE)
haystack <- sample(9500000, 9500000, replace = FALSE)

bench::mark(
  "0.50" = with_load_50(vec_match(needles, haystack)),
  "0.77" = with_load_77(vec_match(needles, haystack)),
  iterations = 20
)
#> # A tibble: 2 × 6
#>   expression      min   median `itr/sec` mem_alloc `gc/sec`
#>   <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
#> 1 0.50          402ms    420ms      2.29     172MB    0.765
#> 2 0.77          499ms    523ms      1.91     108MB    0.478

---

This is typically less useful with functions like vec_unique_loc(), because those only build the dictionary, they don't also use a filled dictionary (which might be near that 77% max capacity) to look up values of some other vector in it, which is typically the painful part.

But it can still be useful when x has a lot of unique values in it, as just building the dictionary will slow down as we approach that 77% load factor.

set.seed(1234)

# 8 mil possible values
dict <- tidyr::expand_grid(
  x = seq_len(200) + 0,
  y = seq_len(200) + 0,
  z = seq_len(200) + 0
)
x <- dplyr::sample_n(dict, size = 1600000, replace = TRUE)

# % unique
vec_unique_count(x) / vec_size(x)
#> [1] 0.9061387

# load factor (max of 77%)
vec_unique_count(x) / 2097152
#> [1] 0.691329

bench::mark(
  "0.50" = with_load_50(vec_unique_loc(x)),
  "0.77" = with_load_77(vec_unique_loc(x)),
  iterations = 10
)
#> # A tibble: 2 × 6
#>   expression      min   median `itr/sec` mem_alloc `gc/sec`
#>   <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
#> 1 0.50          118ms    128ms      7.94    43.6MB     1.98
#> 2 0.77          146ms    151ms      6.52    35.6MB     1.63

It also seems like it doesn't hurt very much in cases where we don't have many unique values (i.e. we don't really approach the max load factor). It does seem to be a little slower.

set.seed(1234)

# Not very many unique values, so nowhere near the max load factor
x <- sample(10, 1600000, replace = TRUE)

bench::mark(
  "0.50" = with_load_50(vec_unique(x)),
  "0.77" = with_load_77(vec_unique(x)),
  iterations = 100
)
#> # A tibble: 2 × 6
#>   expression      min   median `itr/sec` mem_alloc `gc/sec`
#>   <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
#> 1 0.50         6.86ms   9.32ms      107.    22.1MB     8.04
#> 2 0.77         6.75ms    8.3ms      118.    14.1MB     4.91

---

This also helps with the set operations, like vec_set_intersect(), which builds the dictionary on x and looks y up in it. When x has a lot of unique values (which is definitely plausible, considering it is supposed to be a set which doesn't have duplicates), the lower max load factor makes a big difference if the size of x forces us to approach the current 77% max load factor.

And it seems like spending a little bit more memory when we don't actually need it (few uniques) doesn't hurt much.

set.seed(1234)

# set intersection builds hash table on `x`

# 10,000,000 possible values
# 1.6 mil in x, all unique (high load factor for 77% max load case)
# 5.0 mil in y, all unique
x <- sample(1e7, 1600000, replace = FALSE)
y <- sample(1e7, 5000000, replace = FALSE)

bench::mark(
  "0.50" = with_load_50(vec_set_intersect(x, y)),
  "0.77" = with_load_77(vec_set_intersect(x, y)),
  iterations = 100
)
#> # A tibble: 2 × 6
#>   expression      min   median `itr/sec` mem_alloc `gc/sec`
#>   <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
#> 1 0.50          351ms    363ms      2.73    48.8MB    0.683
#> 2 0.77          557ms    624ms      1.59    40.8MB    0.259

# Same idea, but only 10 possible values so load factor is very low.
# Doesn't seem to hurt much.
x <- sample(10, 1600000, replace = TRUE)
y <- sample(10, 5000000, replace = TRUE)

bench::mark(
  "0.50" = with_load_50(vec_set_intersect(x, y)),
  "0.77" = with_load_77(vec_set_intersect(x, y)),
  iterations = 100
)
#> # A tibble: 2 × 6
#>   expression      min   median `itr/sec` mem_alloc `gc/sec`
#>   <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
#> 1 0.50         39.3ms     43ms      22.8    42.7MB     3.11
#> 2 0.77         39.2ms   42.1ms      23.4    34.7MB     3.19

[hadley]

Looks like a nice win!

[lionel-]

Interesting analysis! It seems reasonable to favour speed over space.

[lionel-]

Maybe fetch the load factor from an envvar at load time for easier experimentation?

[lionel-]

oh I guess a "vctrs:::" global option makes more sense

[DavisVaughan]

I think it is reasonable to do this if it comes up again. But for now I don't feel the need to add in the little bit of extra complexity for something that I am hopeful won't need to change again