具有RcppArmadillo的大型SpMat对象

问题描述

我正在尝试学习Rcpp和RcppArmadillo并将其用于稀疏的线性代数例程.

I am trying to learn and use Rcpp and RcppArmadillo for the sparse linear algebra routines.

下面的代码是此处示例的改编: http://gallery.rcpp. org/articles/armadillo-sparse-matrix/

Code below is adaptation of the example here: http://gallery.rcpp.org/articles/armadillo-sparse-matrix/

code <- '
  S4 matx(x);
  IntegerVector Xd = matx.slot("Dim");
  IntegerVector Xi = matx.slot("i");
  IntegerVector Xp = matx.slot("p");
  NumericVector Xx = matx.slot("x");

  arma::sp_mat Xsp(Xd[0], Xd[1]);

  // create space for values, and copy
  arma::access::rw(Xsp.values) = arma::memory::acquire_chunked<double>(Xx.size() + 1);
  arma::arrayops::copy(arma::access::rwp(Xsp.values), 
             Xx.begin(), 
             Xx.size() + 1);

  // create space for row_indices, and copy -- so far in a lame loop
  arma::access::rw(Xsp.row_indices) = arma::memory::acquire_chunked<arma::uword>(Xx.size() + 1);
  for (int j=0; j<Xi.size(); j++)
    arma::access::rwp(Xsp.row_indices)[j] = Xi[j];

  // create space for col_ptrs, and copy -- so far in a lame loop
  arma::access::rw(Xsp.col_ptrs) = arma::memory::acquire_chunked<arma::uword>(Xp.size() + 1);
  for (int j=0; j<Xp.size(); j++)
      arma::access::rwp(Xsp.col_ptrs)[j] = Xp[j];

  // important: set the sentinel as well
  arma::access::rwp(Xsp.col_ptrs)[Xp.size()+1] = std::numeric_limits<arma::uword>::max();

  // set the number of non-zero elements
  arma::access::rw(Xsp.n_nonzero) = Xx.size();

  Rcout << "SpMat Xsp:\\n" << arma::dot(Xsp,Xsp) << std::endl;
'

norm2 <- cxxfunction(signature(x="Matrix"),
                     code,plugin="RcppArmadillo")

当我使用1e4的矢量时，一切正常:

When I use a vector of 1e4, things work fine:

> p <- 10000
> X <- Matrix(rnorm(p),sparse=TRUE)
> norm2(X)
SpMat Xsp:
9997.14
NULL

但是，当我使用长度为1e5的向量时，会产生错误

However, when I use a vector of length 1e5, an error is produced

> p <- 100000
> X <- Matrix(rnorm(p),sparse=TRUE)
> norm2(X)

error: SpMat::init(): requested size is too large

Error: 
>

我似乎无法弄清楚我在做什么错.任何指针将不胜感激.

I cannot seem to figure out what I am doing wrong.Any pointers would be appreciated.

==============更多信息=============

============== more information ==============

问题似乎在于尺寸> = 2 ^ 16 = 65536

The problem seems to be with having dimension >= 2^16=65536

以下作品:

> m <- 1000
> n <- 65535
> nnz <- 10000
> iind <- sample.int(m,nnz,replace=TRUE)
> jind <- sample.int(n,nnz,replace=TRUE)
> xval <- rnorm(nnz)
> X <- sparseMatrix(i=iind,j=jind,x=xval,dims=c(m,n))
> norm2(X)
SpMat Xsp:
10029.8
NULL

以下操作无效:

> m <- 1000
> n <- 65536
> nnz <- 10000
> iind <- sample.int(m,nnz,replace=TRUE)
> jind <- sample.int(n,nnz,replace=TRUE)
> xval <- rnorm(nnz)
> X <- sparseMatrix(i=iind,j=jind,x=xval,dims=c(m,n))
> norm2(X)

error: SpMat::init(): requested size is too large

Error: 
> 

为什么会这样?

推荐答案

您的矩阵似乎很奇怪.通过说

Your matrix seems odd. By saying

 Matrix(rnorm(p),sparse=TRUE)

尽管稀疏，您最终还是得到了一个p x 1的矩阵.如果我只分配10行或列的东西正常工作.

you end up with a p x 1 matrix, albeit sparse. If I just assign 10 rows or colums thingsjust work.

R> p <- 100000
R> X <- Matrix(rnorm(p),nrow=10,sparse=TRUE)
R> dim(X)
[1]    10 10000
R> norm2(X)
SpMat Xsp:
100832
NULL
R> 

所以我认为您只需要一个更好的稀疏矩阵即可使用-转换代码和Armadillo的稀疏矩阵类型就可以了.

So I think you just need a better sparse matrix to work with -- the conversion code, and Armadillo's sparse Matrix type, are fine.

2013年4月30日更新:这实际上是一个Armadillo错误，仅在上游修复. SVN中现在有一个新的RcppArmadillo版本0.3.810.2，应该很快就会迁移到CRAN.您不再需要定义ARMA_64BIT_WORD.

Update on 2013-04-30: This was actually an Armadillo bug, which was just fixed upstream. A new RcppArmadillo verion 0.3.810.2 is now in SVN, and should migrate soon to CRAN shortly. You no longer need to define ARMA_64BIT_WORD.

这篇关于具有RcppArmadillo的大型SpMat对象的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！