block diagonal matrix

Question

Suppose I have a list Z consisting of many matrices and I want to construct a block diagonal matrix from it. ex :

[[1]]
             [,1]        [,2]         [,3]
[1,] 1.002500e+00 0.001930454 1.388794e-11
[2,] 1.930454e-03 1.002500000 1.930454e-03
[3,] 1.388794e-11 0.001930454 1.002500e+00

[[2]]
             [,1]        [,2]         [,3]
[1,] 1.002500e+00 0.001930454 1.388794e-11
[2,] 1.930454e-03 1.002500000 1.930454e-03
[3,] 1.388794e-11 0.001930454 1.002500e+00

I want to create a block diagonal matrix , I am currently using

block = bdiag(z)

However the bdiag command is slow when the number of matrices in the list is large. What is a fast and easy way to construct a block diagonal matrix from the list?

Note my matrix is also symmetric and every matrix in the list has similar dimensions.

I tried with 100,000 3x3 matrices, and it took less than a minute. a list of 200 25x25 matrices, eg, l <- rep(list(matrix(1:625, 25)), 200); x <- Matrix::bdiag(l) finishes almost instantly — rawr
@rawr my application requires optimization of a function which involves bdiag() thus if my bdiag takes 1second then my optimization will take a huge amount of time. i am trying to squeeze the code as much as possible. — raK1
it might be possible to do this faster by indexing into a sparse matrix, but it seems very unlikely that the bdiag() step is really the limiting factor ... ? — Ben Bolker

Stéphane Laurent Stéphane Laurent · Accepted Answer · 2018-07-15T00:42:50

I've attempted to do a faster function, bdiag2.

library(Matrix)

bdiag2 <- function(Ms){
  l <- length(Ms)
  N <- nrow(Ms[[1]])  
  i0 <- rep(1:N, times=N:1)
  s <- rep(seq(0,(l-1)*N,by=N),each=length(i0))
  i <- rep(i0,l) + s
  j0 <- 1:N -> j
  for(k in 1:N){
    j0 <- j0[-1]
    j <- c(j, j0)
  }
  j <- rep(j,l) + s
  idx <- t(upper.tri(Ms[[1]], diag = TRUE))
  x <- unlist(lapply(Ms, "[", idx))
  sparseMatrix(i, j, x = x, symmetric = TRUE)
}

# test & benchmark
N <- 5; l <- 200
Ms <- replicate(l, toeplitz(1:N), simplify = FALSE)

stopifnot(all(bdiag(Ms) == bdiag2(Ms)))

library(microbenchmark)
microbenchmark(
  bdiag = bdiag(Ms),
  bdiag2 = bdiag2(Ms)
)


> microbenchmark(
+   bdiag = bdiag(Ms),
+   bdiag2 = bdiag2(Ms)
+ )
Unit: milliseconds
   expr      min        lq      mean    median        uq       max neval cld
  bdiag 25.68078 27.393898 29.710474 28.566398 30.265242 54.393319   100   b
 bdiag2  1.34185  1.476838  1.693573  1.558724  1.769127  4.482054   100  a

block diagonal matrix

1 Answers