※ 引述《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上来
我或是部分有闲的板友应该可以帮忙看看