11
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ By default, the argument 'replace' of 'sample' is set to FALSE. $\endgroup$ – ocram Mar 8 '12 at 6:29
  • $\begingroup$ Thanks ocram, but that's working within a particular sample. So that ensures that 0,1,2, and 3 won't repeat within a draw (so, I can't draw 0,1,2,2), but I don't know whether that guarantees that the second sample, I can't draw the same sequence of 0123 again. That's what I'm wondering implementation-wise, whether setting the seed has any effect on that repetition. $\endgroup$ – Mittenchops Mar 8 '12 at 20:48
  • $\begingroup$ Yes, this is what I finally understood by reading the answers ;-) $\endgroup$ – ocram Mar 9 '12 at 6:37
  • 1
    $\begingroup$ If limit exceeds 12, you will likely run out of RAM when R attempts to allocate space for seq(1,factorial(limit)). (12! requires about 2 GB, so 13! will need about 25 GB, 14! about 350 GB, etc.) $\endgroup$ – whuber Mar 16 '12 at 19:24
  • 2
    $\begingroup$ There is a fast, compact, and elegant solution for generating random sequences of all permutations of 1:n, provided you can comfortably store n! integers in the range 0:(n!). It combines the inversion table representation of a permutation with factorial base representation of numbers. $\endgroup$ – whuber Mar 16 '12 at 19:31
9
$\begingroup$

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.
$\endgroup$
  • $\begingroup$ This is close, but I notice I get some errors, like 1, 2, and 4, but I think I see what you mean and should be able to work with it. Thanks! > rperm(6,3) $failures [1] 9 $sample [,1] [,2] [,3] [1,] 3 1 3 [2,] 2 2 1 [3,] 1 3 2 [4,] 1 2 2 [5,] 3 3 1 [6,] 2 1 3 $\endgroup$ – Mittenchops Apr 7 '12 at 1:00
3
$\begingroup$

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
$\endgroup$
  • $\begingroup$ Sorry for not explaining the code properly. I'm in a bit of a rush now, but I'm happy to answer any questions you have later. Also, I have no idea about the speed of the above code... $\endgroup$ – MånsT Mar 8 '12 at 8:37
  • 1
    $\begingroup$ I functionalized what you gave me this way: `myperm <- function(limit) { 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)) } return(unique(permutations)) }' It works, but while I can do limit=6, limit=7 makes my computer overheat. =P I think there must still be a way to subsample this... $\endgroup$ – Mittenchops Mar 15 '12 at 5:36
  • $\begingroup$ @Mittenchops, Why do you say we need to use unique for resampling in R without repeating permutations? Thank you. $\endgroup$ – Frank Nov 15 '16 at 4:01
2
$\begingroup$

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
$\endgroup$
  • $\begingroup$ I like this A LOT, and I'm sure it's the right thinking. But my problem has me using a sequence going up to 10. Permn() was significantly slower between factorial(7) and factorial(8), so I think 9 and 10 are going to be prohibitively huge. $\endgroup$ – Mittenchops Mar 8 '12 at 4:57
  • $\begingroup$ @Mittenchops True, but it's still possible you really only need to calculate them once, right? Save them to a file, and then load them when you need them and "sample" from pre-defined list. So you could do the slow calculation of permn(10) or whatever just once. $\endgroup$ – joran Mar 8 '12 at 5:03
  • $\begingroup$ Right, but if I am storing all permutations somewhere, even this breaks down by around factorial(15)---simply too much room to store. That's why I'm wondering whether setting the seed will allow me to sample permutations collectively---and if not, if there is an algorithm for doing that. $\endgroup$ – Mittenchops Mar 8 '12 at 5:16
  • $\begingroup$ @Mittenchops Setting a seed will not influence performance, it just guarantees the same start every time you make a call to PRNG. $\endgroup$ – Roman Luštrik Mar 15 '12 at 8:25
  • 1
    $\begingroup$ @Mitten See the help for set.seed: it describes how to save the state of the RNG and restore it later. $\endgroup$ – whuber Mar 15 '12 at 22:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.