更快地重复初始化 R Matrix::sparseMatrix

huangapple
go评论103阅读模式
英文:

Faster repeated initialization of R Matrix::sparseMatrix

问题

我正在使用R语言处理稀疏矩阵。这些矩阵是方阵，维度接近100万。我需要在代码中重复初始化这些矩阵，而这部分是一个瓶颈。出于各种原因，我使用了i-j-x的表示方法。首先，这些矩阵是块稀疏的，块之间存在冗余结构，使用i-j表示法复制块似乎更容易。此外，我处理的一些矩阵是矩阵的和，仅传递冗余的行索引和列索引比对两个矩阵求和更高效。关于如何构建矩阵，你有什么建议或不建议吗？

编辑：应该像这样：

  1. n_blocks = 1000 # 块的数量
  2. n_within_block = 1000 # 每个块的维度
  3. n_nonnull_within_block = 40000 # 每个块中非零系数的数量
  4. i_within_block = ceiling(n_within_block*runif(n_nonnull_within_block)) # 一个块中的行索引
  5. j_within_block = ceiling(n_within_block*runif(n_nonnull_within_block)) # 一个块中的列索引
  6. x_within_block = rnorm(runif(n_nonnull_within_block))  # 一个块中的系数
  7. # 大多数步骤可以预先计算并重复使用
  8. i = outer( # 行索引
  9.   (ceiling(n_blocks*runif(n_blocks))-1)*n_within_block, 
  10.   i_within_block, 
  11.   "+"
  12. )
  13. j = outer( # 列索引
  14.   (ceiling(n_blocks*runif(n_blocks))-1)*n_within_block, 
  15.   j_within_block, 
  16.   "+"
  17. )
  18. x = rep(x_within_block, n_blocks) # 系数的值
  20. # 优化步骤（希望能提高效率）
  21. t1 = Sys.time()
  22. M = Matrix::sparseMatrix(
  23.   i = i, j = j, x = x 
  24. )
  25. Sys.time()-t1
  26. # 在各个地方使用小的维度！
  27. # Matrix::image(M)
英文:

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

  2. n_blocks = 1000 # number of blocks
  3. n_within_block = 1000 # dimension of each block
  4. n_nonnull_within_block = 40000 # number of non null coefficients within a block
  5. i_within_block = ceiling(n_within_block*runif(n_nonnull_within_block)) # row idx in one block
  6. j_within_block = ceiling(n_within_block*runif(n_nonnull_within_block)) # col idx in one block
  7. x_within_block = rnorm(runif(n_nonnull_within_block))  # coeff in one block
  8. # most of those steps can be pre-computed once and re-used
  9. i = outer( # row idx
  10.   (ceiling(n_blocks*runif(n_blocks))-1)*n_within_block, 
  11.   i_within_block, 
  12.   "+"
  13. )
  14. j = outer( # col idx
  15.   (ceiling(n_blocks*runif(n_blocks))-1)*n_within_block, 
  16.   j_within_block, 
  17.   "+"
  18. )
  19. x =rep(x_within_block, n_blocks) # value of coefficients
  21. # step to optimize (hopefully)
  22. t1 = Sys.time()
  23. M = Matrix::sparseMatrix(
  24.   i = i, j = j, x = x 
  25. )
  26. Sys.time()-t1
  27. # put small dimensions everywhere!
  28. #Matrix::image(M)

答案1

得分: 2

演示了@MikaelJagan的建议，使用整数而不是双精度，并直接构建dgTMatrix：

基准计时：

  1. library(Matrix)
  3. set.seed(1)
  4. n_blocks = 1000 # number of blocks
  5. n_within_block = 1000 # dimension of each block
  6. n_nonnull_within_block = 40000 # number of non null coefficients within a block
  7. i_within_block = ceiling(n_within_block*runif(n_nonnull_within_block)) # row idx in one block
  8. j_within_block = ceiling(n_within_block*runif(n_nonnull_within_block)) # col idx in one block
  9. x_within_block = rnorm(runif(n_nonnull_within_block))  # coeff in one block
  10. # most of those steps can be pre-computed once and re-used
  12. system.time({
  13.   i <- outer( # row idx
  14.     (ceiling(n_blocks*runif(n_blocks))-1)*n_within_block, 
  15.     i_within_block, 
  16.     "+"
  17.   )
  18.   
  19.   j <- outer( # col idx
  20.     (ceiling(n_blocks*runif(n_blocks))-1)*n_within_block, 
  21.     j_within_block, 
  22.     "+"
  23.   )
  24.   
  25.   x <- rep(x_within_block, n_blocks) # value of coefficients
  26.   
  27.   # step to optimize (hopefully)
  28.   M1 <- sparseMatrix(i = i, j = j, x = x, dims = rep(n_blocks*n_within_block, 2))
  29. })
  30. #>    user  system elapsed 
  31. #>    3.28    0.64    3.92

使用整数而不是双精度：

  1. set.seed(1)
  2. n_blocks <- 1000L # number of blocks
  3. n_within_block <- 1000L # dimension of each block
  4. n_nonnull_within_block <- 40000L # number of non null coefficients within a block
  5. i_within_block <- as.integer(n_within_block*runif(n_nonnull_within_block)) # row idx in one block
  6. j_within_block <- as.integer(n_within_block*runif(n_nonnull_within_block)) # col idx in one block
  7. x_within_block <- rnorm(runif(n_nonnull_within_block))  # coeff in one block
  8. # most of those steps can be pre-computed once and re-used
  10. system.time({
  11.   i <- outer( # row idx
  12.     as.integer(n_blocks*runif(n_blocks))*n_within_block, 
  13.     i_within_block, 
  14.     "+"
  15.   )
  16.   j <- outer( # col idx
  17.     as.integer(n_blocks*runif(n_blocks))*n_within_block, 
  18.     j_within_block, 
  19.     "+"
  20.   )
  21.   
  22.   x <- rep(x_within_block, n_blocks) # value of coefficients
  23.   
  24.   M2 <- sparseMatrix(
  25.     i = i, j = j, x = x,
  26.     dims = rep(n_blocks*n_within_block, 2), index1 = FALSE
  27.   )
  28. })
  29. #>    user  system elapsed 
  30. #>    2.41    0.47    2.87

直接构建dgTMatrix：

  1. set.seed(1)
  2. n_blocks <- 1000L # number of blocks
  3. n_within_block <- 1000L # dimension of each block
  4. n_nonnull_within_block <- 40000L # number of non null coefficients within a block
  5. i_within_block <- as.integer(n_within_block*runif(n_nonnull_within_block)) # row idx in one block
  6. j_within_block <- as.integer(n_within_block*runif(n_nonnull_within_block)) # col idx in one block
  7. x_within_block <- rnorm(runif(n_nonnull_within_block))  # coeff in one block
  8. # most of those steps can be pre-computed once and re-used
  10. system.time({
  11.   M3 <- new("dgTMatrix")
  12.   M3@Dim <- rep(n_blocks*n_within_block, 2)
  13.   
  14.   M3@i <- c(outer( # row idx
  15.     as.integer(n_blocks*runif(n_blocks))*n_within_block, 
  16.     i_within_block, 
  17.     "+"
  18.   ))
  19.   
  20.   M3@j <- c(outer( # col idx
  21.     as.integer(n_blocks*runif(n_blocks))*n_within_block, 
  22.     j_within_block, 
  23.     "+"
  24.   ))
  25.   
  26.   M3@x <- rep(x_within_block, n_blocks) # value of coefficients
  27. })
  28. #>    user  system elapsed 
  29. #>    0.66    0.22    0.87

将dgTMatrix转换为CsparseMatrix：

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

请注意，dgTMatrix在内存中比CsparseMatrix占用更多空间：

  1. object.size(M1)
  2. #> 474641504 bytes
  3. object.size(M3)
  4. #> 640001496 bytes
英文:

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

Baseline timing:

  1. library(Matrix)
  3. set.seed(1)
  4. n_blocks = 1000 # number of blocks
  5. n_within_block = 1000 # dimension of each block
  6. n_nonnull_within_block = 40000 # number of non null coefficients within a block
  7. i_within_block = ceiling(n_within_block*runif(n_nonnull_within_block)) # row idx in one block
  8. j_within_block = ceiling(n_within_block*runif(n_nonnull_within_block)) # col idx in one block
  9. x_within_block = rnorm(runif(n_nonnull_within_block))  # coeff in one block
  10. # most of those steps can be pre-computed once and re-used
  12. system.time({
  13.   i &lt;- outer( # row idx
  14.     (ceiling(n_blocks*runif(n_blocks))-1)*n_within_block, 
  15.     i_within_block, 
  16.     &quot;+&quot;
  17.   )
  18.   
  19.   j &lt;- outer( # col idx
  20.     (ceiling(n_blocks*runif(n_blocks))-1)*n_within_block, 
  21.     j_within_block, 
  22.     &quot;+&quot;
  23.   )
  24.   
  25.   x &lt;- rep(x_within_block, n_blocks) # value of coefficients
  26.   
  27.   # step to optimize (hopefully)
  28.   M1 &lt;- sparseMatrix(i = i, j = j, x = x, dims = rep(n_blocks*n_within_block, 2))
  29. })
  30. #&gt;    user  system elapsed 
  31. #&gt;    3.28    0.64    3.92

Using integers instead of doubles:

  1. set.seed(1)
  2. n_blocks &lt;- 1000L # number of blocks
  3. n_within_block &lt;- 1000L # dimension of each block
  4. n_nonnull_within_block &lt;- 40000L # number of non null coefficients within a block
  5. i_within_block &lt;- as.integer(n_within_block*runif(n_nonnull_within_block)) # row idx in one block
  6. j_within_block &lt;- as.integer(n_within_block*runif(n_nonnull_within_block)) # col idx in one block
  7. x_within_block &lt;- rnorm(runif(n_nonnull_within_block))  # coeff in one block
  8. # most of those steps can be pre-computed once and re-used
  10. system.time({
  11.   i &lt;- outer( # row idx
  12.     as.integer(n_blocks*runif(n_blocks))*n_within_block, 
  13.     i_within_block, 
  14.     &quot;+&quot;
  15.   )
  16.   j &lt;- outer( # col idx
  17.     as.integer(n_blocks*runif(n_blocks))*n_within_block, 
  18.     j_within_block, 
  19.     &quot;+&quot;
  20.   )
  21.   
  22.   x &lt;- rep(x_within_block, n_blocks) # value of coefficients
  23.   
  24.   M2 &lt;- sparseMatrix(
  25.     i = i, j = j, x = x,
  26.     dims = rep(n_blocks*n_within_block, 2), index1 = FALSE
  27.   )
  28. })
  29. #&gt;    user  system elapsed 
  30. #&gt;    2.41    0.47    2.87

Building a dgTMatrix directly:

  1. set.seed(1)
  2. n_blocks &lt;- 1000L # number of blocks
  3. n_within_block &lt;- 1000L # dimension of each block
  4. n_nonnull_within_block &lt;- 40000L # number of non null coefficients within a block
  5. i_within_block &lt;- as.integer(n_within_block*runif(n_nonnull_within_block)) # row idx in one block
  6. j_within_block &lt;- as.integer(n_within_block*runif(n_nonnull_within_block)) # col idx in one block
  7. x_within_block &lt;- rnorm(runif(n_nonnull_within_block))  # coeff in one block
  8. # most of those steps can be pre-computed once and re-used
  10. system.time({
  11.   M3 &lt;- new(&quot;dgTMatrix&quot;)
  12.   M3@Dim &lt;- rep(n_blocks*n_within_block, 2)
  13.   
  14.   M3@i &lt;- c(outer( # row idx
  15.     as.integer(n_blocks*runif(n_blocks))*n_within_block, 
  16.     i_within_block, 
  17.     &quot;+&quot;
  18.   ))
  19.   
  20.   M3@j &lt;- c(outer( # col idx
  21.     as.integer(n_blocks*runif(n_blocks))*n_within_block, 
  22.     j_within_block, 
  23.     &quot;+&quot;
  24.   ))
  25.   
  26.   M3@x &lt;- rep(x_within_block, n_blocks) # value of coefficients
  27. })
  28. #&gt;    user  system elapsed 
  29. #&gt;    0.66    0.22    0.87

Converting a dgTMatrix to a CsparseMatrix:

  1. system.time(M4 &lt;- as(M3, &quot;CsparseMatrix&quot;))
  2. #&gt;    user  system elapsed 
  3. #&gt;    1.47    0.20    1.68
  4. identical(M1, M2)
  5. #&gt; [1] TRUE
  6. identical(M1, M4)
  7. #&gt; [1] TRUE

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

  1. object.size(M1)
  2. #&gt; 474641504 bytes
  3. object.size(M3)
  4. #&gt; 640001496 bytes

huangapple
  • 本文由 发表于 2023年8月9日 17:42:30
  • 转载请务必保留本文链接：https://go.coder-hub.com/76866469.html
  • r
  • sparse-matrix
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定