在R中更新稀疏矩阵的单个元素非常缓慢 - 如何更快地进行更新?

3
在使用稀疏矩阵的脚本中,profvis分析表明,稀疏矩阵元素的更新是该过程中最慢的步骤,慢1个数量级。我需要了解是否可以更好(尤其是更快),希望有人能指导我去哪里找或提供建议。这里是一些R代码,可以重现我脚本的“关键”部分:
require(Matrix)
m <- new("dgCMatrix", i = c(0L, 1L, 2L, 6L, 8L, 0L, 1L, 2L, 5L, 6L, 
7L, 0L, 1L, 2L, 7L, 3L, 4L, 3L, 4L, 5L, 6L, 1L, 4L, 5L, 6L, 0L, 
1L, 4L, 5L, 6L, 8L, 10L, 1L, 2L, 7L, 0L, 6L, 8L, 9L, 10L, 6L, 
9L, 10L), p = c(0L, 5L, 11L, 15L, 17L, 21L, 25L, 32L, 35L, 38L, 
40L, 43L), Dim = c(11L, 11L), Dimnames = list(c("1", "2", "3", 
"4", "5", "6", "7", "8", "9", "10", "11"), c("1", "2", "3", "4", 
"5", "6", "7", "8", "9", "10", "11")), x = c(2, 1, 1, 1, 1, 1, 
3, 2, 1, 1, 1, 1, 2, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 2, 1, 1, 
1, 2, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2), factors = list())

system.time(for (i in 1:10000) m[7,c(1,5,6,7)] <- c(0,0,1,0))

在我的笔记本上,这大约需要7秒钟。 [顺便说一句,显然我没有重复执行同一操作10000次;更新的行和列每次都会改变,但确实会发生大量的更新操作。我之前这么做是为了模拟在真实脚本中执行的操作,并获得可以与更快的解决方案进行比较的可测时间。]
有什么想法/建议吗?
附言 我过去也遇到过类似的问题,但是是不同的情况;由于我的活动历史只能追溯到几个月前,所以找不回来了。
编辑 好的,我找到了如何检索所有旧帖子的方法,并发现我这里描述的问题并未被涉及。
编辑2 - 跟进/pseudospin的建议
require(Matrix)
require(data.table)

m <- new("dgCMatrix", i = c(0L, 1L, 2L, 6L, 8L, 0L, 1L, 2L, 5L, 6L, 
                            7L, 0L, 1L, 2L, 7L, 3L, 4L, 3L, 4L, 5L, 6L, 1L, 4L, 5L, 6L, 0L, 
                            1L, 4L, 5L, 6L, 8L, 10L, 1L, 2L, 7L, 0L, 6L, 8L, 9L, 10L, 6L, 
                            9L, 10L), p = c(0L, 5L, 11L, 15L, 17L, 21L, 25L, 32L, 35L, 38L, 
                                            40L, 43L), Dim = c(11L, 11L), Dimnames = list(c("1", "2", "3", 
                                                                                            "4", "5", "6", "7", "8", "9", "10", "11"), c("1", "2", "3", "4", 
                                                                                                                                         "5", "6", "7", "8", "9", "10", "11")), x = c(2, 1, 1, 1, 1, 1, 
                                                                                                                                                                                      3, 2, 1, 1, 1, 1, 2, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 2, 1, 1, 
                                                                                                                                                                                      1, 2, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2), factors = list())

ms <- summary(m)
ms <- ms[order(ms$i,ms$j),]
msdt <- data.table(ms)

time_1 <- system.time(for (i in 1:5000) m[7,c(1,5,7,9)] <- c(0,0,1,0))
cat("\ntime_1 =", time_1)

time_2 <- system.time(for (i in 1:5000) ms[(ms$i == 7) & (ms$j %in% c(1,5,7,9)),"x"] <- c(0,0,1,0))
cat("\ntime_2 =", time_2)

time_3 <- system.time(for (i in 1:5000) msdt[(i == 7) & (j %in% c(1,5,7,9)),"x" := c(0,0,1,0)])
cat("\ntime_3 =", time_3)

这给出了:

time_1 = 2.86 0 2.86 NA NA
time_2 = 0.23 0 0.24 NA NA
time_3 = 1.2 0.02 1.22 NA NA

也许这个例子有点误导人,因为通常我将有更高的i和 j 的最大值,所以从data.table 子集取数据可能比从 data.frame 子集取数据更有效率。
需要使用我的真实数据进行测试...

编辑3 - 使用GKi建议的密集矩阵方法尝试真实数据测试

真实数据(太大无法粘贴在此):m 是一个稀疏的 5828 x 5828 的矩阵;其中 302986/33965584 = 0.9% 是填充的(因此是稀疏的)。它占用4.4 MB。对应的密集矩阵 dm = as.matrix(m) 占用272.5MB。

测试更新方法 sparseMatrix (1), data.frame (2), data.table (3) 和密集矩阵 (4) 显示如下:

time_1 = 10.25 3.19 13.72 NA NA
time_2 = 41.32 10.94 52.52 NA NA
time_3 = 35.64 7.44 43.34 NA NA
time_4 = 0.05 0.03 0.08 NA NA

因此,与GKi的结果一致,密集矩阵法远远是最快的方法,但代价是巨大的存储空间。
另一方面,最初使用的模拟数据对于方法给出了非常不同的结果,实际数据表明该方法在这4种方法中排名第二。

不幸的是,这看起来像一个进退两难的局面:为了进行快速编辑,我需要使用密集矩阵,但密集矩阵需要太多的内存,因此我需要使用稀疏矩阵,但是编辑速度慢 :(

也许我需要重新考虑pseudospin的原始建议,并使用矩阵每行的稀疏向量。为此,我需要找出如何通过间接引用(字符串)引用存储的R对象。

我猜这个程序运行“慢”是因为稀疏矩阵的存储方式 - 它必须查找每个元素以确定是否存在,然后再进行更改,而传统矩阵只需使用从索引计算出的偏移量找到内存进行更改。我认为你的问题可能需要一个更全面的解决方案来使它更快。 - pseudospin
我明白了,谢谢。我在想,是否用整个矩阵行替换而不是替换行中的一些元素会更好呢? - user6376297
2
这里有一个关于稀疏矩阵存储方式的不错介绍。基于此,直接编辑稀疏矩阵并不高效 - 每次编辑子集元素(无论是一个元素还是一行)都可能会导致完全复制。 - pseudospin
也许如果我使用稀疏矩阵汇总数据框(i,j,x)并对其进行子集和编辑,效果会更好?也许可以考虑使用数据表以提高速度? - user6376297
不确定使用data.table实现可编辑的稀疏矩阵是否比Matrix包更好。可能需要考虑实际项目来找出最佳方法,但是猜测,将行存储为单独的稀疏向量可能是高效的方法? - pseudospin
显示剩余5条评论
1个回答

2

不同方法的比较。请注意,DataFrameDataTableFastMatch这些方法仅适用于在稀疏矩阵中覆盖现有值时使用,而不是在插入新值时使用。

library(microbenchmark)
microbenchmark(list = fun, control=list(order="block"))
#Unit: microseconds
#         expr     min       lq      mean   median       uq      max neval   cld
# SparseMatrix 618.307 623.1795 639.38102 627.2525 643.3895 1021.301   100     e
#  DenseMatrix   1.178   1.2210   1.36957   1.2635   1.3435    7.060   100 a    
#         Slam 259.703 264.3945 270.57151 265.9780 268.0745  426.610   100   c  
#        Spray 422.129 427.1310 463.21071 434.0705 440.6025 2880.787   100    d 
#    DataFrame  37.031  37.7910  38.98143  38.1660  38.6255   73.283   100  b   
#   DataFrameB  16.928  17.4480  17.85553  17.6910  18.0155   28.859   100 ab   
#   DataFrameC  21.007  21.7170  22.68689  21.9600  22.3735   38.175   100 ab   
#    DataTable 283.409 288.5710 299.43498 292.8395 301.1255  500.484   100   c  
#    FastMatch  40.138  40.7885  42.33623  41.2165  41.6575   82.274   100  b   
#         List   2.703   2.7900   3.28163   2.8535   2.9770   13.375   100 a    
#  Environment   2.157   2.2055   2.29915   2.2575   2.3340    4.211   100 a    

library(bench)
mark(exprs = fun, check = FALSE)
## A tibble: 11 x 13
#   expression        min   median `itr/sec` mem_alloc `gc/sec` n_itr  n_gc
#   <bch:expr>   <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl> <int> <dbl>
# 1 SparseMatrix 624.43µs 651.98µs     1487.    1.65KB    10.4    712     5
# 2 DenseMatrix    1.61µs   1.86µs   490003.        0B    49.0   9999     1
# 3 Slam         261.46µs 271.96µs     3606.     7.2KB     8.24  1750     4
# 4 Spray        424.08µs 439.86µs     2206.    8.16KB    10.4   1063     5
# 5 DataFrame     37.26µs  40.03µs    24378.    2.08KB     9.75  9996     4
# 6 DataFrameB    18.25µs  19.72µs    49447.     1.7KB     9.89  9998     2
# 7 DataFrameC    22.72µs  24.82µs    39142.      840B    11.7   9997     3
# 8 DataTable    288.24µs 300.25µs     3252.   18.34KB     8.24  1579     4
# 9 FastMatch     41.05µs  43.73µs    22292.    2.46KB    11.2   9995     5
#10 List            3.4µs   3.71µs   257225.        0B    25.7   9999     1
#11 Environment    2.82µs   3.11µs   306445.        0B     0    10000     0

方法:

fun <- alist(
  SparseMatrix = m[7,c(1,5,6,7)] <- c(0,0,1,0)
, DenseMatrix = dm[7,c(1,5,6,7)] <- c(0,0,1,0)
, Slam = slm[7,c(1,5,6,7)] <- c(0,0,1,0)
, Spray = spm[7,c(1,5,6,7)] <- c(0,0,1,0)
, DataFrame = ms[(ms$j == 7) & (ms$i %in% c(1,5,6,7)),"x"] <- c(0,0,1,0)
, DataFrameB = ms$x[(ms$j == 7) & (ms$i %in% c(1,5,6,7))] <- c(0,0,1,0)
, DataFrameC = {i <- which(ms$j == 7); ms$x[i[ms$i[i] %in% c(1,5,6,7)]] <- c(0,0,1,0)}
, DataTable = msdt[(j == 7) & (i %in% c(1,5,6,7)),"x" := c(0,0,1,0)]
, FastMatch = mf[mf$j %fin% 7 & (mf$i %fin% c(1,5,6,7)),"x"] <- c(0,0,1,0)
, List = ml[["7"]][c("1","5","6","7")] <- c(0,0,1,0)
, Environment = me[["7"]][c("1","5","6","7")] <- c(0,0,1,0)
)

数据:

library(Matrix)
m <- new("dgCMatrix", i = c(0L, 1L, 2L, 6L, 8L, 0L, 1L, 2L, 5L, 6L, 
7L, 0L, 1L, 2L, 7L, 3L, 4L, 3L, 4L, 5L, 6L, 1L, 4L, 5L, 6L, 0L, 
1L, 4L, 5L, 6L, 8L, 10L, 1L, 2L, 7L, 0L, 6L, 8L, 9L, 10L, 6L, 
9L, 10L), p = c(0L, 5L, 11L, 15L, 17L, 21L, 25L, 32L, 35L, 38L, 
40L, 43L), Dim = c(11L, 11L), Dimnames = list(c("1", "2", "3", 
"4", "5", "6", "7", "8", "9", "10", "11"), c("1", "2", "3", "4", 
"5", "6", "7", "8", "9", "10", "11")), x = c(2, 1, 1, 1, 1, 1, 
3, 2, 1, 1, 1, 1, 2, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 2, 1, 1, 
1, 2, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2), factors = list())

dm <- as.matrix(m) #Dense Matrix

library(slam)
slm <- as.simple_sparse_array(dm)

library(spray)
spm <- as.spray(dm)

ms <- summary(m)
ms <- ms[order(ms$i,ms$j),]

library(data.table)
msdt <- data.table(ms)

library(fastmatch)
mf <- ms

ml <- split(setNames(ms$x, ms$j), ms$i)

me <- list2env(ml, hash = TRUE)

非常感谢您进行的全面分析!我希望我可以使用一个密集矩阵...但很遗憾,由于我所处理的问题规模较大,一个密集矩阵会导致R崩溃(内存不足)。为了给您一个想法,我是通过将(MxN)矩阵与其转置的叉积来获得它的,其中M可以是10000、50000甚至数百万。如果你试图创建一个(100万×100万)的密集矩阵,你可以想象会发生什么。它的大部分元素都是零,因此使用稀疏矩阵是不可避免的。我将用我的真实数据尝试您的基准测试,并报告结果。 - user6376297
我已经添加了 data.frame 的变种版本,还添加了使用 listenvironment 的版本。 - GKi
(゚o゚) - 我不知道这个 split 方法!这几乎和密集矩阵一样好!不确定它是如何工作的,但如果可以实现,它比我迄今见过的任何东西都要好。 - user6376297
还要看看环境,因为它通过哈希查找其元素,而不像列表那样进行线性搜索。 - GKi

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接