※ 引述《chu1216 (chu)》之銘言：
: 請問一下
: 我想要找非零的index的矩陣,
: 因此我用which(XXX != 0, arr.ind = T),
: 但因為矩陣的size非常大, 跑起來花很長時間,
: 請問有類似且效率比較好的的function嗎?
: 感謝!!
我試了一下RcppArmadillo 做了一個小試驗
從結果來看，可能C++可以幫上一點忙，但是要看你sparse的比例
程式碼:
library(Rcpp)
sourceCpp(code = "
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
// [[Rcpp::export]]
arma::uvec find_nonzero_index(const arma::mat & x) {
return find(x != 0) + 1;
}")
genSparseMat <- function(matSize, sparseProp) {
x <- matrix(runif(prod(matSize)), matSize[1])
x[x >= sparseProp] <- 0
return(x)
}
set.seed(100)
matSp30 <- genSparseMat(c(5000, 5000), 0.3)
matSp20 <- genSparseMat(c(5000, 5000), 0.2)
matSp10 <- genSparseMat(c(5000, 5000), 0.1)
matSp05 <- genSparseMat(c(5000, 5000), 0.05)
matSp01 <- genSparseMat(c(5000, 5000), 0.01)
library(microbenchmark)
microbenchmark(
RWhichSp30 = which(matSp30 != 0),
RcppSp30 = find_nonzero_index(matSp30),
RWhichSp20 = which(matSp20 != 0),
RcppSp20 = find_nonzero_index(matSp20),
RWhichSp10 = which(matSp10 != 0),
RcppSp10 = find_nonzero_index(matSp10),
RWhichSp05 = which(matSp05 != 0),
RcppSp05 = find_nonzero_index(matSp05),
RWhichSp01 = which(matSp01 != 0),
RcppSp01 = find_nonzero_index(matSp01),
times = 20L
)
結果:
Unit: milliseconds
expr min lq mean median uq max neval
RWhichSp30 167.1455 182.53765 195.99211 189.95580 207.34010 254.4168 20
RcppSp30 120.1717 125.33950 134.36280 128.67385 134.05210 216.7681 20
RWhichSp20 141.7320 148.71530 164.81187 159.17785 173.29780 224.4695 20
RcppSp20 92.8394 94.95830 99.97492 96.96545 102.32200 122.8803 20
RWhichSp10 118.1888 127.71755 138.10948 136.37605 150.66480 162.9312 20
RcppSp10 53.6464 55.00425 59.48229 59.18855 63.43540 68.5329 20
RWhichSp05 106.1757 111.50920 127.11906 117.38275 133.84795 231.4787 20
RcppSp05 35.1256 36.11950 38.13195 38.11350 39.82205 41.8294 20
RWhichSp01 95.0594 102.33750 113.32983 107.19800 124.01145 150.2782 20
RcppSp01 19.9087 20.64855 21.85202 21.43630 22.75725 26.0881 20
可能的方向:
1. 計算中，產生matrix的時候，是否就可以直接是SparseMatrix，而不用轉換
ex: 用Matrix::sparse.model.matrix 而非使用model.matrix
2. 整個計算移到RcppArmadillo的架構，透過C++加速迴圈 (R迴圈很慢)
3. 重新設計演算法，避免取大量非零元素的index
4. 寫一個OpenMP的C++函數幫你用多個cores跑
因為你沒有提供其他細節，所以我只能提供這樣的方向建議
如果你有一個最小可重現的程式的話，歡迎PO上來
我或是部分有閒的板友應該可以幫忙看看