How to resample in R without repeating permutations?

Question

In R, if I set.seed(), and then use the sample function to randomize a list, can I guarantee I won't generate the same permutation?

ie...

set.seed(25)
limit <- 3
myindex <- seq(0,limit)
for (x in seq(1,factorial(limit))) {
    permutations <- sample(myindex)
    print(permutations)
}

This produces

[1] 1 2 0 3
[1] 0 2 1 3
[1] 0 3 2 1
[1] 3 1 2 0
[1] 2 3 0 1
[1] 0 1 3 2

will all permutations printed be unique permutations? Or is there some chance, based on the way this is implemented, that I could get some repeats?

I want to be able to do this without repeats, guaranteed. How would I do that?

(I also want to avoid having to use a function like permn(), which has a very mechanistic method for generating all permutations---it doesn't look random.)

Also, sidenote---it looks like this problem is O((n!)!), if I'm not mistaken.

Answer 1

The question has many valid interpretations. The comments--especially the one indicating permutations of 15 or more elements are needed (15! = 1307674368000 is getting big)--suggest that what is wanted is a relatively small random sample, without replacement, of all n! = n*(n-1)*(n-2)*...*2*1 permutations of 1:n. If this is true, there exist (somewhat) efficient solutions.

The following function, rperm, accepts two arguments n (the size of the permutations to sample) and m (the number of permutations of size n to draw). If m approaches or exceeds n!, the function will take a long time and return many NA values: it is intended for use when n is relatively big (say, 8 or more) and m is much smaller than n!. It works by caching a string representation of the permutations found so far and then generating new permutations (randomly) until a new one is found. It exploits R's associative list-indexing ability to search the list of previously-found permutations quickly.

rperm <- function(m, size=2) { # Obtain m unique permutations of 1:size

    # Function to obtain a new permutation.
    newperm <- function() {
        count <- 0                # Protects against infinite loops
        repeat {
            # Generate a permutation and check against previous ones.
            p <- sample(1:size)
            hash.p <- paste(p, collapse="")
            if (is.null(cache[[hash.p]])) break

            # Prepare to try again.
            count <- count+1
            if (count > 1000) {   # 1000 is arbitrary; adjust to taste
                p <- NA           # NA indicates a new permutation wasn't found
                hash.p <- ""
                break
            }
        }
        cache[[hash.p]] <<- TRUE  # Update the list of permutations found
        p                         # Return this (new) permutation
    }

    # Obtain m unique permutations.
    cache <- list()
    replicate(m, newperm())  
} # Returns a `size` by `m` matrix; each column is a permutation of 1:size.

The nature of replicate is to return the permutations as column vectors; e.g., the following reproduces an example in the original question, transposed:

> set.seed(17)
> rperm(6, size=4)
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    2    4    4    3    4
[2,]    3    4    1    3    1    2
[3,]    4    1    3    2    2    3
[4,]    2    3    2    1    4    1

Timings are excellent for small to moderate values of m, up to about 10,000, but degrade for larger problems. For example, a sample of m = 10,000 permutations of n = 1000 elements (a matrix of 10 million values) was obtained in 10 seconds; a sample of m = 20,000 permutations of n = 20 elements required 11 seconds, even though the output (a matrix of 400,000 entries) was much smaller; and computing sample of m = 100,000 permutations of n = 20 elements was aborted after 260 seconds (I didn't have the patience to wait for completion). This scaling problem appears to be related to scaling inefficiencies in R's associative addressing. One can work around it by generating samples in groups of, say, 1000 or so, then combining those samples into a large sample and removing duplicates. R experts might be able to suggest more efficient solutions or better workarounds.

Edit

We can achieve near linear asymptotic performance by breaking the cache into a hierarchy of two caches, so that R never has to search through a large list. Conceptually (although not as implemented), create an array indexed by the first $k$ elements of a permutation. Entries in this array are lists of all permutations sharing those first $k$ elements. To check whether a permutation has been seen, use its first $k$ elements to find its entry in the cache and then search for that permutation within that entry. We can choose $k$ to balance the expected sizes of all the lists. The actual implementation does not use a $k$-fold array, which would be hard to program in sufficient generality, but instead uses another list.

Here are some elapsed times in seconds for a range of permutation sizes and numbers of distinct permutations requested:

 Number Size=10 Size=15 Size=1000 size=10000 size=100000
     10    0.00    0.00      0.02       0.08        1.03
    100    0.01    0.01      0.07       0.64        8.36
   1000    0.08    0.09      0.68       6.38
  10000    0.83    0.87      7.04      65.74
 100000   11.77   10.51     69.33
1000000  195.5   125.5

(The apparently anomalous speedup from size=10 to size=15 is because the first level of the cache is larger for size=15, reducing the average number of entries in the second-level lists, thereby speeding up R's associative search. At some cost in RAM, execution could be made faster by increasing the upper-level cache size. Just increasing k.head by 1 (which multiplies the upper-level size by 10) sped up rperm(100000, size=10) from 11.77 seconds to 8.72 seconds, for instance. Making the upper-level cache 10 times bigger yet achieved no appreciable gain, clocking at 8.51 seconds.)

Except for the case of 1,000,000 unique permutations of 10 elements (a substantial portion of all 10! = about 3.63 million such permutations), practically no collisions were ever detected. In this exceptional case, there were 169,301 collisions, but no complete failures (one million unique permutations were in fact obtained).

Note that with large permutation sizes (greater than 20 or so), the chance of obtaining two identical permutations even in a sample as large as 1,000,000,000 is vanishingly small. Thus, this solution is applicable primarily in situations where (a) large numbers of unique permutations of (b) between $n=5$ and $n=15$ or so elements are to be generated but even so, (c) substantially fewer than all $n!$ permutations are needed.

Working code follows.

rperm <- function(m, size=2) { # Obtain m unique permutations of 1:size
    max.failures <- 10

    # Function to index into the upper-level cache.
    prefix <- function(p, k) {    # p is a permutation, k is the prefix size
        sum((p[1:k] - 1) * (size ^ ((1:k)-1))) + 1
    } # Returns a value from 1 through size^k

    # Function to obtain a new permutation.
    newperm <- function() {
        # References cache, k.head, and failures in parent context.
        # Modifies cache and failures.        

        count <- 0                # Protects against infinite loops
        repeat {
            # Generate a permutation and check against previous ones.
            p <- sample(1:size)
            k <- prefix(p, k.head)
            ip <- cache[[k]]
            hash.p <- paste(tail(p,-k.head), collapse="")
            if (is.null(ip[[hash.p]])) break

            # Prepare to try again.
            n.failures <<- n.failures + 1
            count <- count+1
            if (count > max.failures) {  
                p <- NA           # NA indicates a new permutation wasn't found
                hash.p <- ""
                break
            }
        }
        if (count <= max.failures) {
            ip[[hash.p]] <- TRUE      # Update the list of permutations found
            cache[[k]] <<- ip
        }
        p                         # Return this (new) permutation
    }

    # Initialize the cache.
    k.head <- min(size-1, max(1, floor(log(m / log(m)) / log(size))))
    cache <- as.list(1:(size^k.head))
    for (i in 1:(size^k.head)) cache[[i]] <- list()

    # Count failures (for benchmarking and error checking).
    n.failures <- 0

    # Obtain (up to) m unique permutations.
    s <- replicate(m, newperm())
    s[is.na(s)] <- NULL
    list(failures=n.failures, sample=matrix(unlist(s), ncol=size))
} # Returns an m by size matrix; each row is a permutation of 1:size.

Answer 2

Using unique in the right way ought to do the trick:

set.seed(2)
limit <- 3
myindex <- seq(0,limit)

endDim<-factorial(limit)
permutations<-sample(myindex)

while(is.null(dim(unique(permutations))) || dim(unique(permutations))[1]!=endDim) {
    permutations <- rbind(permutations,sample(myindex))
}
# Resulting permutations:
unique(permutations)

# Compare to
set.seed(2)
permutations<-sample(myindex)
for(i in 1:endDim)
{
permutations<-rbind(permutations,sample(myindex))
}
permutations
# which contains the same permutation twice

Answer 3

I"m going to side step your first question a bit, and suggest that if your are dealing with relatively short vectors, you could simply generate all the permutations using permn and them randomly order those using sample:

x <- combinat:::permn(1:3)
> x[sample(factorial(3),factorial(3),replace = FALSE)]
[[1]]
[1] 1 2 3

[[2]]
[1] 3 2 1

[[3]]
[1] 3 1 2

[[4]]
[1] 2 1 3

[[5]]
[1] 2 3 1

[[6]]
[1] 1 3 2

Answer 4

The package RcppAlgos (I am the author) has some functions specifically for this type of problem. Observe:

library(RcppAlgos)
permuteSample(v = 0:10, m = 5, seed = 123, n = 3)
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]   10    2    5    6    3
#> [2,]    0    6   10    3    5
#> [3,]    5   10    3    2    9

## Or use global RNG to get the same results
set.seed(123)
permuteSample(v = 0:10, m = 5, n = 3)
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]   10    2    5    6    3
#> [2,]    0    6   10    3    5
#> [3,]    5   10    3    2    9

They are guaranteed to sample without replacement:

permuteCount(3)
#> [1] 6

permuteSample(3, n = 7)
#> Error: n exceeds the maximum number of possible results

If we request the maximal number of samples, it will be equivalent to generating all results and shuffling:

ps <- permuteSample(3, n = 6, seed = 100)

## Shuffled
ps
#>      [,1] [,2] [,3]
#> [1,]    1    3    2
#> [2,]    2    1    3
#> [3,]    1    2    3
#> [4,]    3    2    1
#> [5,]    3    1    2
#> [6,]    2    3    1

## Same as getting them all
identical(
    permuteGeneral(3),
    ps[do.call(order, as.data.frame(ps)), ]
)
#> [1] TRUE

Flexibility

These functions are quite general as well. For example, if we want to sample a vector where repetition within the vector is allowed set repetition = TRUE. This is not to be confused with sampling with replacement (see the section: Sampling With Replacement below):

permuteSample(v = 0:10, m = 5, repetition = TRUE, seed = 42, n = 3)
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]    4    2    1    5   10
#> [2,]    3    7    9    8    7
#> [3,]    6    8    8    8    5

## Go beyond the vector length
permuteSample(v = 3, m = 10, repetition = TRUE, seed = 42, n = 3)
#>      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#> [1,]    3    3    1    3    2    3    3    3    1     2
#> [2,]    2    3    1    2    2    3    3    3    3     3
#> [3,]    1    2    2    1    1    1    3    3    2     3

What about specific multiplicity for each element… easy! Use the freqs parameter:

## With v = 3 and freqs = 2:4, this means:
##
##    1 can occur a maximum of 2 times
##    2 can occur a maximum of 3 times
##    3 can occur a maximum of 4 times
##
permuteSample(v = 3, m = 6, freqs = 2:4, seed = 1234, n = 3)
#>      [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,]    3    1    3    3    2    3
#> [2,]    3    2    3    2    2    1
#> [3,]    3    3    3    2    3    2

It works on all atomic types as well:

## When m is NULL, the length defaults to the length of v
set.seed(28)
permuteSample(
    v = rnorm(3) + rnorm(3) * 1i,
    repetition = TRUE,
    n = 5
)
#>                      [,1]                   [,2]                   [,3]
#> [1,] -1.331167+0.5313963i -1.33116707+0.5313963i -1.90215722-1.8199917i
#> [2,] -1.902157-1.8199917i -1.33116707+0.5313963i -0.06429479+0.1626697i
#> [3,] -1.331167+0.5313963i -0.06429479+0.1626697i -1.33116707+0.5313963i
#> [4,] -1.902157-1.8199917i -1.90215722-1.8199917i -1.33116707+0.5313963i
#> [5,] -1.902157-1.8199917i -1.90215722-1.8199917i -1.90215722-1.8199917i

Sampling With Replacement

The default behavior of the sampling functions in RcppAlgos is to sample without replacement. If you need sampling with replacement, we can make use of the sampleVec parameter. This parameter takes a vector and interprets them as indices representing the lexicographical permutations to return. For example:

permuteSample(3, sampleVec = c(1, 3, 6))
#>      [,1] [,2] [,3]
#> [1,]    1    2    3
#> [2,]    2    1    3
#> [3,]    3    2    1

## Same as
permuteGeneral(3)[c(1, 3, 6), ]
#>      [,1] [,2] [,3]
#> [1,]    1    2    3
#> [2,]    2    1    3
#> [3,]    3    2    1

This means we can generate our own random sample of indices and pass it to permuteSample:

set.seed(123)
idx <- sample(permuteCount(3), 5, replace = TRUE)
idx
#> [1] 3 6 3 2 2

permuteSample(3, sampleVec = idx)
#>      [,1] [,2] [,3]
#> [1,]    2    1    3
#> [2,]    3    2    1
#> [3,]    2    1    3
#> [4,]    1    3    2
#> [5,]    1    3    2

## See rownames corresponding to the lexicographical result 
permuteSample(3, sampleVec = idx, namedSample = TRUE)
#>   [,1] [,2] [,3]
#> 3    2    1    3
#> 6    3    2    1
#> 3    2    1    3
#> 2    1    3    2
#> 2    1    3    2

Efficiency

The functions are written in C++ with efficiency in mind. Take for example the largest case @whuber tests:

system.time(size_10 <- permuteSample(10, n = 1e6, seed = 42))
#>    user  system elapsed 
#>   0.269   0.006   0.275

## Guaranteed no repeated permutations!
nrow(unique(size_10))
#> [1] 1000000

system.time(size_15 <- permuteSample(15, n = 1e6, seed = 42))
#>    user  system elapsed 
#>   0.395   0.008   0.403

## Guaranteed no repeated permutations!
nrow(unique(size_15))
#> [1] 1000000

Why not try even larger cases?!?!?! Using nThreads, we can achieve even greater efficiency:

## Single threaded
system.time(permuteSample(15, n = 1e7, seed = 321))
#>    user  system elapsed 
#>   4.490   0.108   4.599

## Using 4 threads
system.time(permuteSample(15, n = 1e7, seed = 321, nThreads = 4))
#>    user  system elapsed 
#>   4.794   0.083   2.106

Sampling Other Combinatorial Objects

There are also sampling functions for other combinatorial objects:

1. permutations: permuteSample
2. combinations: comboSample
3. integer partitions: partitionsSample
4. integer compositions: compositionsSample
5. partitions of a vector into groups comboGroupsSample