What is the fastest way to create a block diagonal matrix other than bdiag function in the matrix package which is slow for large datasets I have provided a reproducible example below library(Matrix) cy11=function(a,b,H){ d=outer(a,b,`-`);I=outer(a,b,`==`) H[1]^2*d^2/H[2]^2+I*H[3]^2 } x0=c(1,1,0.05) z=list() for (i in 1:100) { a=rnorm(50,0,10) b=rnorm(50,0,10) z[[i]]=cy11(a,a,x0) } bdiag(z) [[alternative HTML version deleted]]