计算带状矩阵的 colCumsums 的更快替代方法

Question

我是 R 和 stats 的新手。在我目前工作的领域中，我需要以独特的方式计算累积列总和。

最初提供一个宽度为 b 且行数为 n 的方带矩阵。例如对于 n = 8 和 b = 3

0 1 2 7 0 0 0 0
0 0 3 6 7 0 0 0
0 0 0 3 1 7 0 0
0 0 0 0 4 4 7 0
0 0 0 0 0 5 8 7
0 0 0 0 0 0 1 8
0 0 0 0 0 0 0 4
0 0 0 0 0 0 0 0   

然后对矩阵进行变换，得到一个以对角线为列的 n x b 矩阵。就像给定的例子一样，

1 2 7  
3 6 7 
3 1 7 
4 4 7 
5 8 7 
1 8 0
4 0 0
0 0 0

我目前正在使用以下函数来执行此操作。

     packedband <- function(x, n, b) {
      mat <- sapply(0:(b-1), function(i)
         diag(x[-(n:(n-i)), -(1:(1+i))])[1:n] )
      mat[is.na(mat)] <- 0
      return(mat)
      }

然后应用 matrixStats 包中的 colCumsums 函数来获得所需的输出矩阵。对于给定的示例，

1    2     7
4    8    14
7    9    21
11   13   28
16   21   35
17   29   35
21   29   35
21   29   35

我正在寻找的是更快地计算这些操作，因为在给定的域中，列（或行）的数量可以> 10^5。从最终目标开始，可能可以删除计算打包带函数的步骤是获取累积列总和。 提前致谢。

Answer 1

在处理了稀疏矩阵之后，我认为for 循环在这里可能会很好用。

尝试原始数据

d = as.matrix(read.table(text="0 1 2 7 0 0 0 0
0 0 3 6 7 0 0 0
0 0 0 3 1 7 0 0
0 0 0 0 4 4 7 0
0 0 0 0 0 5 8 7
0 0 0 0 0 0 1 8
0 0 0 0 0 0 0 4
0 0 0 0 0 0 0 0 "))

colnames(d) <- NULL

功能

packedband <- function(x, b=3) {
      n = nrow(d)
      mat <- sapply(0:(b-1), function(i)
                  diag(x[-(n:(n-i)), -(1:(1+i))])[1:n] )
      mat[is.na(mat)] <- 0
      matrixStats::colCumsums(mat)
   }

forloop <- function(d, b=3){
     n = nrow(d)
     m = matrix(0, n, b)
      for(i in 1:b) {
        ro = 1:(n-i)
        co = (1+i):n
        vec = `length<-`(d[cbind(ro, co)], n)
        vec[is.na(vec)] <- 0
        m[ , i] = cumsum(vec)
      }
     m
   }

# create initial sparse matrix just to omit time to convert
# as if its faster it may be worth storing your band matrices in sparse format
library(Matrix)
m <- as(d, "TsparseMatrix") 

spm <- function(m, b=3){
x = sparseMatrix(i = m@i+1,
                 j = m@j - m@i,
                 x = m@x,
                 dims = c(nrow(m),b))
matrixStats::colCumsums(as.matrix(x))
}

all.equal(forloop(d), packedband(d))
all.equal(spm(m), packedband(d))

尝试更大的数据

d = matrix(0, 5e3, 5e3)
d[(col(d) - row(d)) == 1] <- 1
d[(col(d) - row(d)) == 2] <- 1
d[ (col(d) - row(d)) == 3] <- 1

m <- as(d, "TsparseMatrix") 

all.equal(forloop(d), packedband(d))
all.equal(spm(m), packedband(d))

microbenchmark::microbenchmark(packedband(d), forloop(d), spm(m), times=50)
# Unit: microseconds
#           expr         min          lq        mean      median          uq         max neval cld
#  packedband(d) 1348240.520 1724714.293 1740634.707 1733305.192 1763377.869 1960353.263    50   b
#     forloop(d)     720.344     973.658    1054.461    1026.807    1174.731    1565.912    50  a 
#         spm(m)    2145.875    2437.321    2586.503    2480.133    2749.019    3766.051    50  a