Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Taking the transpose of square blocks in a rectangular matrix r

Suppose I have two square matrices (actually many more) that are bound together:

mat = matrix(1:18,nrow=3,ncol=6)

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

I want to take the transpose of each (3x3) matrix and keep them glued side by side, so the result is:

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

I do not want to do this manually because it is MANY matrices cbound together, not just 2.

I would like a solution that avoids looping or apply (which is just a wrapper for a loop). I need the efficient solution because this will have to run tens of thousands of times.

like image 262
wolfsatthedoor Avatar asked May 30 '16 20:05

wolfsatthedoor


2 Answers

One way is to use matrix indexing

matrix(t(m), nrow=nrow(m))[, c(matrix(1:ncol(m), nrow(m), byrow=T)) ]

This takes the transposed matrix and rearanges the columns in the desired order.

m <- matrix(1:18,nrow=3,ncol=6)
matrix(t(m), nrow=nrow(m))
#      [,1] [,2] [,3] [,4] [,5] [,6]
# [1,]    1   10    2   11    3   12
# [2,]    4   13    5   14    6   15
# [3,]    7   16    8   17    9   18

So we want the 1st, 3rd, and 5th columns, and 2, 4, and 6th columns together. One way is to index these with

c(matrix(1:ncol(m), nrow(m), byrow=T))
#[1] 1 3 5 2 4 6

As an alternative, you could use

idx <- rep(1:ncol(m), each=nrow(m), length=ncol(m)) ;
do.call(cbind, split.data.frame(t(m), idx))

Try on a new matrix

(m <- matrix(1:50, nrow=5))
#      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
# [1,]    1    6   11   16   21   26   31   36   41    46
# [2,]    2    7   12   17   22   27   32   37   42    47
# [3,]    3    8   13   18   23   28   33   38   43    48
# [4,]    4    9   14   19   24   29   34   39   44    49
# [5,]    5   10   15   20   25   30   35   40   45    50

matrix(t(m), nrow=nrow(m))[, c(matrix(1:ncol(m), nrow(m), byrow=T)) ]
#      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
# [1,]    1    2    3    4    5   26   27   28   29    30
# [2,]    6    7    8    9   10   31   32   33   34    35
# [3,]   11   12   13   14   15   36   37   38   39    40
# [4,]   16   17   18   19   20   41   42   43   44    45
# [5,]   21   22   23   24   25   46   47   48   49    50
like image 70
user20650 Avatar answered Oct 24 '22 20:10

user20650


This might do it:

mat = matrix(1:18,nrow=3,ncol=6)
mat

output <- lapply(seq(3, ncol(mat), 3), function(i) { t(mat[, c((i - 2):i)]) } )
output

do.call(cbind, output)

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

Was curious and timed the two approaches. The matrix approach used by user20650 is much faster than the lapply approach I used:

library(microbenchmark)

mat = matrix(1:1600, nrow=4, byrow = FALSE)

lapply.function <- function(x) {

     step1 <- lapply(seq(nrow(mat), ncol(mat), nrow(mat)), function(i) {
                     t(mat[, c((i - (nrow(mat) - 1) ):i)])
              } )

     l.output <- do.call(cbind, step1)
     return(l.output)
}

lapply.output <- lapply.function(mat)

matrix.function <- function(x) { 
     m.output <- matrix(t(mat), nrow=nrow(mat))[, c(matrix(1:ncol(mat), nrow(mat), byrow=TRUE)) ] 
}

matrix.output <- matrix.function(mat)

identical(lapply.function(mat), matrix.function(mat))

microbenchmark(lapply.function(mat), matrix.function(mat), times = 1000)

#Unit: microseconds
#                 expr     min      lq      mean  median      uq      max neval
# lapply.function(mat) 735.602 776.652 824.44917 791.443 809.856 2260.834  1000
# matrix.function(mat)  32.298  35.619  37.75495  36.826  37.732   78.481  1000
like image 41
Mark Miller Avatar answered Oct 24 '22 21:10

Mark Miller