Return all matrices of a given dimension with a certain number of ones and remaining zeros

Viewed 169

Consider the following simplified example, with all possible 2 x 2 matrices with one 1 and the remaining 0s.

library(arrangements)

# define function
generate_matrices <- function(nrow, ncol, ones_count) {
  
  vectors <- permutations(c(
    rep(1, ones_count),
    rep(0, nrow * ncol - ones_count)
  ))
  
  # remove redundancies
  vectors <- vectors[!duplicated(vectors),]
  
  # list of matrices
  out <- list()
  
  for (i in 1:ncol(vectors)) {
    out[[i]] <- matrix(
      data = vectors[,i],
      nrow = nrow,
      ncol = ncol,
      byrow = TRUE
    )
  }
  return(out)
}

Run function to generate all 2 by 2 matrices with one 1:

generate_matrices(nrow = 2, ncol = 2, ones_count = 1)

[[1]]
     [,1] [,2]
[1,]    1    0
[2,]    0    0

[[2]]
     [,1] [,2]
[1,]    0    1
[2,]    0    0

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

[[4]]
     [,1] [,2]
[1,]    0    0
[2,]    0    1

When I extend this to a matrix with 5 rows, 4 columns and 4 ones, it errors:

generate_matrices(nrow = 5, ncol = 4, ones_count = 4)
# Error in permutations(c(rep(1, ones_count), rep(0, nrow * ncol - ones_count))) :
# too many results

My guess is that the lines

vectors <- permutations(c(
    rep(1, ones_count),
    rep(0, nrow * ncol - ones_count)
  ))

takes too long to run and/or there is not enough memory on my laptop to store these. Is there a more efficient way to implement this?

It is worth noting that I would like to eventually extend this to the 6 x 5 case with 4 ones, and 8 x 5 case with 8 ones.

4 Answers

You can take combination of indices on which is 1:

m <- 2
n <- 2
k <- 2

createMatrix <- function(m, n, indices){
  
  x <- matrix(0, m, n)
  x[indices] <- 1
  
  x
}

lapply(
  combn(seq_len(m*n), k, simplify = FALSE), 
  function(x) createMatrix(m, n, x)
)

where m is number of rows, n number of columns and k number of ones.

You could use the function developed here to get all possible matrices of dim 5x4, and then filter by the number of 1s using sum.

f = function(nrow, ncol) lapply(asplit(do.call(expand.grid, rep(list(0:1), nrow * ncol)), 1), matrix, nrow, ncol)
list = f(5,4)
list[lapply(list, sum) == 4]

Using partitions::multiset and convert result to array of appropriate dimensions seems more efficient.

f1 = function(nr, nc, n1){
  m = unclass(multiset(c(rep(0, nr*nc - n1), rep(1, n1))))
  `dim<-`(m, c(nr, nc, ncol(m)))
}

f1(nr = 2, nc = 2, n1 = 1)
# , , 1,
# 
# [,1] [,2]
# [1,]    0    0
# [2,]    0    1
# 
# , , 2
# 
# [,1] [,2]
# [1,]    0    1
# [2,]    0    0
# 
# , , 3
# 
# [,1] [,2]
# [1,]    0    0
# [2,]    1    0
# 
# , , 4
# 
# [,1] [,2]
# [1,]    1    0
# [2,]    0    0

If needed, the array can easily be converted to a list:

asplit(a, MARGIN = 3) 

Benchmark

On larger data, multiset is considerably faster than the answer using combn (which needs to call matrix and [<- for every combination). Here timings on 8*5 matrices with 6 ones, which results in 3 838 380 matrices:

f2 = function(m, n, k){lapply(
  combn(seq_len(m*n), k, simplify = FALSE), 
  function(x) createMatrix(m, n, x))}

microbenchmark(
  f1(nr = 8, nc = 5, n1 = 6),
  f2(m = 8, n = 5, k = 6),
  times = 10L)
    
# Unit: milliseconds
#                         expr        min         lq       mean     median        uq       max neval
#  f1(nr = 8, nc = 5, n1 = 6)    582.5020   680.5886   916.1864   802.1724 1531.4137  3132.456    10
#      f2(m = 8, n = 5, k = 6) 20539.4030 22039.6975 24097.4683 24022.0033 1166.9455  2544.132    10

dim(f1(nr = 8, nc = 5, n1 = 6))
# [1] 8 5 3838380

length(f2(m = 8, n = 5, k = 6))
# [1] 3838380

With the input above, the code by Maël unfortunately errored on my PC ("cannot allocate vector of size..."), possibly due to the "expand.grid explosion".

Here is a one-liner using the package RcppAlgos (I am the author):

library(RcppAlgos)
nr = 2
nc = 2
n1 = 1

permuteGeneral(0:1, freqs = c(nr * nc - n1, n1),
               FUN = function(x) matrix(x, nc))
[[1]]
     [,1] [,2]
[1,]    0    0
[2,]    0    1

[[2]]
     [,1] [,2]
[1,]    0    1
[2,]    0    0

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

[[4]]
     [,1] [,2]
[1,]    1    0
[2,]    0    0

This package also offers very flexible iterators via permuteIter. For example:

iter <- permuteIter(0:1, freqs = c(nr * nc - n1, n1),
                    FUN = function(x) matrix(x, nc))
iter$nextIter()
     [,1] [,2]
[1,]    0    0
[2,]    0    1

iter$back() ## Get the last one (or the first one via front())
     [,1] [,2]
[1,]    1    0
[2,]    0    0

iter[[3]] ## Random access via [[
     [,1] [,2]
[1,]    0    0
[2,]    1    0

iter[[c(1, 4)]]
[[1]]
     [,1] [,2]
[1,]    0    0
[2,]    0    1

[[2]]
     [,1] [,2]
[1,]    1    0
[2,]    0    0

iter$startOver() ## Reset the iterator
iter$nextNIter(3) ## Get the next n
[[1]]
     [,1] [,2]
[1,]    0    0
[2,]    0    1

[[2]]
     [,1] [,2]
[1,]    0    1
[2,]    0    0

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

iter$prevIter() ## iterate backwards
     [,1] [,2]
[1,]    0    1
[2,]    0    0

The package is written in C++ and generally very efficient, however using the FUN parameter is purely for convenience.

If you are looking for raw speed, you are better off using the approach offered up by @Henrik. That is, generate all permutations as a matrix and convert it to an array. Note the dimension is in a different order. This is so, because partitions::multiset generates a column-major object while RcppAlgos::permuteGeneral generates a row-major matrix. The resulting arrays are isomorphic.

fasterApproach <- function(nr, nc, n1) {
    arr <- permuteGeneral(0:1, freqs = c(nr * nc - n1, n1))
    dim(arr) <- c(nrow(arr), nr, nc)
    arr
}

dim(f1(6, 5, 5))
[1]      6      5 142506

dim(fasterApproach(6, 5, 5))
[1] 142506      6      5

microbenchmark(f1(6, 5, 5), fasterApproach(6, 5, 5))
Unit: milliseconds
                    expr       min       lq     mean   median       uq      max neval cld
             f1(6, 5, 5) 14.240662 22.09967 34.19180 24.15314 28.84547 125.2753   100   b
 fasterApproach(6, 5, 5)  9.006603 10.15762 20.41521 15.87181 18.20326 115.5324   100  a
Related