3
votes

I would like to create a block diagonal matrix with the diagonal blocks being repeated a certain number of times, with the off-diagonal blocks being all zero matrices. For example, suppose we start with a matrix with the following:

> diag.matrix

       [,1] [,2] [,3] [,4] [,5]
  [1,]  1.0  0.5  0.5  0.5  0.5
  [2,]  0.5  1.0  0.5  0.5  0.5
  [3,]  0.5  0.5  1.0  0.5  0.5
  [4,]  0.5  0.5  0.5  1.0  0.5
  [5,]  0.5  0.5  0.5  0.5  1.0

I would like this matrix to be the diagonal block matrix so that in the end I have something like:

       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]  1.0  0.5  0.5  0.5  0.5  0.0  0.0  0.0  0.0   0.0
 [2,]  0.5  1.0  0.5  0.5  0.5  0.0  0.0  0.0  0.0   0.0
 [3,]  0.5  0.5  1.0  0.5  0.5  0.0  0.0  0.0  0.0   0.0
 [4,]  0.5  0.5  0.5  1.0  0.5  0.0  0.0  0.0  0.0   0.0
 [5,]  0.5  0.5  0.5  0.5  1.0  0.0  0.0  0.0  0.0   0.0
 [6,]  0.0  0.0  0.0  0.0  0.0  1.0  0.5  0.5  0.5   0.5
 [7,]  0.0  0.0  0.0  0.0  0.0  0.5  1.0  0.5  0.5   0.5
 [8,]  0.0  0.0  0.0  0.0  0.0  0.5  0.5  1.0  0.5   0.5
 [9,]  0.0  0.0  0.0  0.0  0.0  0.5  0.5  0.5  1.0   0.5
[10,]  0.0  0.0  0.0  0.0  0.0  0.5  0.5  0.5  0.5   1.0

Here, we have the same block matrix above repeated twice in the block diagonals. If I wanted to do this efficiently an arbitrary number of times, is there a way to do it? thanks!

1

1 Answers

10
votes

1) kronecker If M is your matrix and k is the number of times you want it repeated then:

kronecker(diag(k), M)

For example,

M <- matrix(0.5, 5, 5) + diag(0.5, 5)
k <- 2
kronecker(diag(k), M)

giving:

      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]  1.0  0.5  0.5  0.5  0.5  0.0  0.0  0.0  0.0   0.0
 [2,]  0.5  1.0  0.5  0.5  0.5  0.0  0.0  0.0  0.0   0.0
 [3,]  0.5  0.5  1.0  0.5  0.5  0.0  0.0  0.0  0.0   0.0
 [4,]  0.5  0.5  0.5  1.0  0.5  0.0  0.0  0.0  0.0   0.0
 [5,]  0.5  0.5  0.5  0.5  1.0  0.0  0.0  0.0  0.0   0.0
 [6,]  0.0  0.0  0.0  0.0  0.0  1.0  0.5  0.5  0.5   0.5
 [7,]  0.0  0.0  0.0  0.0  0.0  0.5  1.0  0.5  0.5   0.5
 [8,]  0.0  0.0  0.0  0.0  0.0  0.5  0.5  1.0  0.5   0.5
 [9,]  0.0  0.0  0.0  0.0  0.0  0.5  0.5  0.5  1.0   0.5
[10,]  0.0  0.0  0.0  0.0  0.0  0.5  0.5  0.5  0.5   1.0

1a) %x% The last line of code can alternately be written as:

diag(k) %x% M

2) Matrix::bdiag Another possibility if you want to save space is to create a sparse matrix of class "dgMCatrix". It does not store the zero values. See ?bdiag :

library(Matrix)

bdiag(replicate(k, M, simplify = FALSE))

giving:

10 x 10 sparse Matrix of class "dgCMatrix"

 [1,] 1.0 0.5 0.5 0.5 0.5 .   .   .   .   .  
 [2,] 0.5 1.0 0.5 0.5 0.5 .   .   .   .   .  
 [3,] 0.5 0.5 1.0 0.5 0.5 .   .   .   .   .  
 [4,] 0.5 0.5 0.5 1.0 0.5 .   .   .   .   .  
 [5,] 0.5 0.5 0.5 0.5 1.0 .   .   .   .   .  
 [6,] .   .   .   .   .   1.0 0.5 0.5 0.5 0.5
 [7,] .   .   .   .   .   0.5 1.0 0.5 0.5 0.5
 [8,] .   .   .   .   .   0.5 0.5 1.0 0.5 0.5
 [9,] .   .   .   .   .   0.5 0.5 0.5 1.0 0.5
[10,] .   .   .   .   .   0.5 0.5 0.5 0.5 1.0

2b) Diagonal or to create a sparse matrix of class "dgTMatrix" :

Diagonal(k) %x% M