Faster repeated initialization of R Matrix::sparseMatrix

I am working with sparse matrices in R. They are square, with dimensions neighboring 1 million. I have to initialize those matrices repeatedly in my code, and this part is a bottleneck. I am using the i-j-x formulation, for various reasons. First, the matrices are block-sparse, with redundant structures across blocks, and it seems infinitely easier to replicate the blocks using the i-j formulation. Also, some of the matrices I work with are sums of matrices, and it is much more efficient just to pass redundant is and js than to sum two matrices. Any tips about how I should or should not build my matrices ?

Edit: should look like that

n_blocks = 1000 # number of blocks
n_within_block = 1000 # dimension of each block
n_nonnull_within_block = 40000 # number of non null coefficients within a block
i_within_block = ceiling(n_within_block*runif(n_nonnull_within_block)) # row idx in one block
j_within_block = ceiling(n_within_block*runif(n_nonnull_within_block)) # col idx in one block
x_within_block = rnorm(runif(n_nonnull_within_block))  # coeff in one block
# most of those steps can be pre-computed once and re-used
i = outer( # row idx
  (ceiling(n_blocks*runif(n_blocks))-1)*n_within_block, 
  i_within_block, 
  "+"
)
j = outer( # col idx
  (ceiling(n_blocks*runif(n_blocks))-1)*n_within_block, 
  j_within_block, 
  "+"
)
x =rep(x_within_block, n_blocks) # value of coefficients

# step to optimize (hopefully)
t1 = Sys.time()
M = Matrix::sparseMatrix(
  i = i, j = j, x = x 
)
Sys.time()-t1
# put small dimensions everywhere!
#Matrix::image(M)

Demonstrating @MikaelJagan's suggestions of using integers instead of doubles and building a dgTMatrix directly:

Baseline timing:

library(Matrix)

set.seed(1)
n_blocks = 1000 # number of blocks
n_within_block = 1000 # dimension of each block
n_nonnull_within_block = 40000 # number of non null coefficients within a block
i_within_block = ceiling(n_within_block*runif(n_nonnull_within_block)) # row idx in one block
j_within_block = ceiling(n_within_block*runif(n_nonnull_within_block)) # col idx in one block
x_within_block = rnorm(runif(n_nonnull_within_block))  # coeff in one block
# most of those steps can be pre-computed once and re-used

system.time({
  i <- outer( # row idx
    (ceiling(n_blocks*runif(n_blocks))-1)*n_within_block, 
    i_within_block, 
    "+"
  )
  
  j <- outer( # col idx
    (ceiling(n_blocks*runif(n_blocks))-1)*n_within_block, 
    j_within_block, 
    "+"
  )
  
  x <- rep(x_within_block, n_blocks) # value of coefficients
  
  # step to optimize (hopefully)
  M1 <- sparseMatrix(i = i, j = j, x = x, dims = rep(n_blocks*n_within_block, 2))
})
#>    user  system elapsed 
#>    3.28    0.64    3.92

Using integers instead of doubles:

set.seed(1)
n_blocks <- 1000L # number of blocks
n_within_block <- 1000L # dimension of each block
n_nonnull_within_block <- 40000L # number of non null coefficients within a block
i_within_block <- as.integer(n_within_block*runif(n_nonnull_within_block)) # row idx in one block
j_within_block <- as.integer(n_within_block*runif(n_nonnull_within_block)) # col idx in one block
x_within_block <- rnorm(runif(n_nonnull_within_block))  # coeff in one block
# most of those steps can be pre-computed once and re-used

system.time({
  i <- outer( # row idx
    as.integer(n_blocks*runif(n_blocks))*n_within_block, 
    i_within_block, 
    "+"
  )
  j <- outer( # col idx
    as.integer(n_blocks*runif(n_blocks))*n_within_block, 
    j_within_block, 
    "+"
  )
  
  x <- rep(x_within_block, n_blocks) # value of coefficients
  
  M2 <- sparseMatrix(
    i = i, j = j, x = x,
    dims = rep(n_blocks*n_within_block, 2), index1 = FALSE
  )
})
#>    user  system elapsed 
#>    2.41    0.47    2.87

Building a dgTMatrix directly:

set.seed(1)
n_blocks <- 1000L # number of blocks
n_within_block <- 1000L # dimension of each block
n_nonnull_within_block <- 40000L # number of non null coefficients within a block
i_within_block <- as.integer(n_within_block*runif(n_nonnull_within_block)) # row idx in one block
j_within_block <- as.integer(n_within_block*runif(n_nonnull_within_block)) # col idx in one block
x_within_block <- rnorm(runif(n_nonnull_within_block))  # coeff in one block
# most of those steps can be pre-computed once and re-used

system.time({
  M3 <- new("dgTMatrix")
  M3@Dim <- rep(n_blocks*n_within_block, 2)
  
  M3@i <- c(outer( # row idx
    as.integer(n_blocks*runif(n_blocks))*n_within_block, 
    i_within_block, 
    "+"
  ))
  
  M3@j <- c(outer( # col idx
    as.integer(n_blocks*runif(n_blocks))*n_within_block, 
    j_within_block, 
    "+"
  ))
  
  M3@x <- rep(x_within_block, n_blocks) # value of coefficients
})
#>    user  system elapsed 
#>    0.66    0.22    0.87

Converting a dgTMatrix to a CsparseMatrix:

system.time(M4 <- as(M3, "CsparseMatrix"))
#>    user  system elapsed 
#>    1.47    0.20    1.68
identical(M1, M2)
#> [1] TRUE
identical(M1, M4)
#> [1] TRUE

Note that the dgTMatrix is larger in memory than the CsparseMatrix:

object.size(M1)
#> 474641504 bytes
object.size(M3)
#> 640001496 bytes