快速替代計算帶矩陣的colCumsums

我是R和stats的新手。在當前工作的域中，我需要以獨特的方式計算累積列總和。快速替代計算帶矩陣的colCumsums

最初寬度b和行數n的正方形帶矩陣是對於n = 8和b provided.For例= 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

則矩陣是在這樣一種方式，anxb待轉化用對角線爲列矩陣是obtained.Like對於給定的例子中，

我目前使用下面的函數來執行此操作。

 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功能packageto獲得期望的輸出matrix.For給出的示例，

我所尋找的是這些操作的，因爲在給定一個更快的計算域，列（或行）的數量可以> 10^5。可能的是，由於最終目標是獲得累積列總和，因此可以去除計算打包帶功能的步驟。在此先感謝。

來源

2017-04-02 Mkty

我不明白，應該是什麼樣的輸出試試。如果你只是想使用colSums（）函數？還有函數'abandSparse（n，m = n，k，diagonals，symmetric = FALSE，giveCsparse = TRUE）'。 – Mislav

@mkty'colCumsums'不是一個基本的R函數。請提及該函數所在的包的名稱。 @Mislav與'abandSparse'相同。 – lmo

@lmo我修改了這個問題。請看看。 – Mkty

在使用稀疏矩陣搞亂之後，我認爲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 = [email protected]+1, 
       j = [email protected] - [email protected], 
       x = [email protected], 
       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

來源

2017-04-03 12:59:16 user20650

這是一個答案的寶石。非常感謝這個解決方案！ – Mkty

快速替代計算帶矩陣的colCumsums

回答

相關問題