4

I have a vector given and want to transform it into a certain block matrix. Consider this simple example:

k <- c(1,2,3)
a <- rep(apply(expand.grid(k, k), 1, prod), each=3)
a
[1] 1 1 1 2 2 2 3 3 3 2 2 2 4 4 4 6 6 6 3 3 3 6 6 6 9 9 9

This vector should be aligned in a block matrix of the form:

rbind(
cbind(diag(a[1:3]), diag(a[4:6]), diag(a[7:9])),
cbind(diag(a[10:12]), diag(a[13:15]), diag(a[16:18]) ),
cbind(diag(a[19:21]), diag(a[22:24]), diag(a[25:27]) ) 
)

       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
 [1,]    1    0    0    2    0    0    3    0    0
 [2,]    0    1    0    0    2    0    0    3    0
 [3,]    0    0    1    0    0    2    0    0    3
 [4,]    2    0    0    4    0    0    6    0    0
 [5,]    0    2    0    0    4    0    0    6    0
 [6,]    0    0    2    0    0    4    0    0    6
 [7,]    3    0    0    6    0    0    9    0    0
 [8,]    0    3    0    0    6    0    0    9    0
 [9,]    0    0    3    0    0    6    0    0    9

Now this is of course a small and simple example and I would like to do this for much larger vectors/matrices. Thus my question: Is there a general way to align a vector in a block matrix of certain form (without looping)?

Sotos
  • 51,121
  • 6
  • 32
  • 66
user436994
  • 601
  • 5
  • 15
  • Maybe relevant post: https://stackoverflow.com/questions/17495841/block-diagonal-binding-of-matrices – zx8754 Oct 20 '17 at 09:33

2 Answers2

6

Instead of manually doing the split, we can use %/%

k <- 3
lst <- split(a, (seq_along(a)-1)%/%k + 1)
do.call(rbind, lapply(split(lst, (seq_along(lst)-1) %/% k + 1), 
       function(x) do.call(cbind,  lapply(x, function(y) diag(y)))))
#       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
# [1,]    1    0    0    2    0    0    3    0    0
# [2,]    0    1    0    0    2    0    0    3    0
# [3,]    0    0    1    0    0    2    0    0    3
# [4,]    2    0    0    4    0    0    6    0    0
# [5,]    0    2    0    0    4    0    0    6    0
# [6,]    0    0    2    0    0    4    0    0    6
# [7,]    3    0    0    6    0    0    9    0    0
# [8,]    0    3    0    0    6    0    0    9    0
# [9,]    0    0    3    0    0    6    0    0    9
Axeman
  • 32,068
  • 8
  • 81
  • 94
akrun
  • 874,273
  • 37
  • 540
  • 662
3

An alternative using the Kronecker product on a slightly different vector is as follows.

# create initial vector
aNew <- rep(1:3, 3) * rep(1:3, each=3)
aNew
[1] 1 2 3 2 4 6 3 6 9

Note that aNew is the unique values of the vector a in the same order, that is, it is equivalent to unique(a). Convert aNew into a 3X3 matrix and then perform the Kronecker product against it and the 3X3 identity matrix.

matrix(aNew, 3, 3) %x% diag(3)
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
 [1,]    1    0    0    2    0    0    3    0    0
 [2,]    0    1    0    0    2    0    0    3    0
 [3,]    0    0    1    0    0    2    0    0    3
 [4,]    2    0    0    4    0    0    6    0    0
 [5,]    0    2    0    0    4    0    0    6    0
 [6,]    0    0    2    0    0    4    0    0    6
 [7,]    3    0    0    6    0    0    9    0    0
 [8,]    0    3    0    0    6    0    0    9    0
 [9,]    0    0    3    0    0    6    0    0    9
lmo
  • 37,904
  • 9
  • 56
  • 69