10

I want to split matrices of size k x l into blocks of size n x n considering an ofset o (Like Mathematica's Partition function does).

For example, given a matrix A like

A <- matrix(seq(1:16), nrow = 4, ncol = 4)

     [,1] [,2] [,3] [,4]
[1,]    1    5    9   13
[2,]    2    6   10   14
[3,]    3    7   11   15
[4,]    4    8   12   16

and block size = 3, offset = 1, I want as output the four submatrices that I'd get from

A[1:3, 1:3]
A[1:3, 2:4]
A[2:4, 1:3]
A[2:4, 2:4]

If offset were equal to 2 or 3, the output for this example should be only the submatrix that I get from

A[1:3, 1:3]

How can I vectorize this?

989
  • 12,579
  • 5
  • 31
  • 53
andandandand
  • 21,946
  • 60
  • 170
  • 271
  • You could index negatively, though this may be case-specific: `apply(combn(c(-1, -1, -4, -4), 2), 2, function(x){list(A[x[1], x[2]])})` – alistaire Jun 19 '16 at 03:49
  • What should we have if `offset == 2`? – Psidom Jun 19 '16 at 04:06
  • @Psidom, an offset of 2 rows and 2 columns. In the 4 x 4 matrix of the example above offsets 2 and 3 both give the submatrix `A[1:3, 1:3]` – andandandand Jun 19 '16 at 04:09
  • 1
    This might be similar - http://stackoverflow.com/questions/24299171/function-to-split-a-matrix-into-sub-matrices-in-r/24300097 - particularly if the simplification-to-array parts are removed – thelatemail Jun 19 '16 at 11:02

2 Answers2

3

There might be a more elegant way. Here is how I'd do it by writing a myPartition function which simulates the mathematica Partition function. Firstly use Map to construct possible index along the row and column axis where we use seq to take offset into consideration, and then use cross2 from purrr to construct a list of all possible combinations of the subset index. Finally use lapply to subset the matrix and return a list of subset matrix;

The testing results on offset 1, 2 and 3 are as follows which seems to behave as expected:

library(purrr)
ind <- function(k, n, o) Map(`:`, seq(1, k-n+1, by = o), seq(n, k, by = o))

# this is a little helper function that generates subset index according to dimension of the 
# matrix, the first sequence construct the starting point of the subset index with an interval 
# of o which is the offset while the second sequence construct the ending point of the subset index
# use Map to construct vector from start to end which in OP's case will be 1:3 and 2:4. 

myPartition <- function(mat, n, o) {
    lapply(cross2(ind(nrow(mat),n,o), ind(ncol(mat),n,o)), function(i) mat[i[[1]], i[[2]]])
}

# This is basically an lapply. we use cross2 to construct combinations of all subset index
# which will be 1:3 and 1:3, 1:3 and 2:4, 2:4 and 1:3 and 2:4 and 2:4 in OP's case. Use lapply
# to loop through the index and subset.

# Testing case for offset = 1
myPartition(A, 3, 1)

# [[1]]
#      [,1] [,2] [,3]
# [1,]    1    5    9
# [2,]    2    6   10
# [3,]    3    7   11

# [[2]]
#      [,1] [,2] [,3]
# [1,]    2    6   10
# [2,]    3    7   11
# [3,]    4    8   12

# [[3]]
#      [,1] [,2] [,3]
# [1,]    5    9   13
# [2,]    6   10   14
# [3,]    7   11   15

# [[4]]
#      [,1] [,2] [,3]
# [1,]    6   10   14
# [2,]    7   11   15
# [3,]    8   12   16

# Testing case for offset = 2
myPartition(A, 3, 2)
# [[1]]
#      [,1] [,2] [,3]
# [1,]    1    5    9
# [2,]    2    6   10
# [3,]    3    7   11

# Testing case for offset = 3
myPartition(A, 3, 3)
# [[1]]
#      [,1] [,2] [,3]
# [1,]    1    5    9
# [2,]    2    6   10
# [3,]    3    7   11
Psidom
  • 209,562
  • 33
  • 339
  • 356
2

How about this using base R, the idea is to generate all possible windows (i.e. winds) of size n*n while taking into account the offset. Then print all possible permutations of winds's elements in matrix A (i.e. perms). It works for any A of size k*l.

A <- matrix(seq(1:16), nrow = 4, ncol = 4)
c <- ncol(A); r <- nrow(A)
offset <- 1; size <- 3
sq <- seq(1, max(r,c), offset)
winds <- t(sapply(sq, function(x) c(x,(x+size-1))))
winds <- winds[winds[,2]<=max(r, c),] # check the range
if (is.vector(winds)) dim(winds) <- c(1,2) # vector to matrix
perms <- expand.grid(list(1:nrow(winds), 1:nrow(winds)))
out=apply(perms, 1, function(x) {
   a11 <- winds[x[1],1];a12 <- winds[x[1],2];a21 <- winds[x[2],1];a22 <- winds[x[2],2]
   if (ifelse(r<c, a12<=r, a22<=c)) { # check the range
       cat("A[", a11, ":", a12, ", ", a21, ":", a22, "]", sep="", "\n")
       print(A[a11:a12, a21:a22])
   }
})

# A[1:3, 1:3]
     # [,1] [,2] [,3]
# [1,]    1    5    9
# [2,]    2    6   10
# [3,]    3    7   11
# A[2:4, 1:3]
     # [,1] [,2] [,3]
# [1,]    2    6   10
# [2,]    3    7   11
# [3,]    4    8   12
# A[1:3, 2:4]
     # [,1] [,2] [,3]
# [1,]    5    9   13
# [2,]    6   10   14
# [3,]    7   11   15
# A[2:4, 2:4]
     # [,1] [,2] [,3]
# [1,]    6   10   14
# [2,]    7   11   15
# [3,]    8   12   16

For size=3 and offset=2 or offset=3:

# A[1:3, 1:3]
     # [,1] [,2] [,3]
# [1,]    1    5    9
# [2,]    2    6   10
# [3,]    3    7   11

For offset=2 and size=2:

# A[1:2, 1:2]
     # [,1] [,2]
# [1,]    1    5
# [2,]    2    6
# A[3:4, 1:2]
     # [,1] [,2]
# [1,]    3    7
# [2,]    4    8
# A[1:2, 3:4]
     # [,1] [,2]
# [1,]    9   13
# [2,]   10   14
# A[3:4, 3:4]
     # [,1] [,2]
# [1,]   11   15
# [2,]   12   16
989
  • 12,579
  • 5
  • 31
  • 53