带有 RcppArmadillo 的大型 SpMat 对象

Question

我正在尝试学习和使用 Rcpp 和 RcppArmadillo 来处理稀疏线性代数例程。

以下代码改编自此处的示例：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 的向量时，一切正常：

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

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

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

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

Error: 
>

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

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

问题似乎在于尺寸 >= 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: 
> 

为什么会这样？

Answer 1

你的矩阵看起来很奇怪。通过说

 Matrix(rnorm(p),sparse=TRUE)

你最终得到一个 p x 1 矩阵，尽管它很稀疏。如果我只分配 10 行或列的东西 只是工作。

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> 

所以我认为你只需要一个更好的稀疏矩阵来使用——转换代码和犰狳的稀疏矩阵类型都很好。

2013-04-30 更新：这实际上是一个 Armadillo 错误，刚刚在上游修复。新的 RcppArmadillo 版本 0.3.810.2 现在在 SVN 中，应该很快就会迁移到 CRAN。您不再需要定义ARMA_64BIT_WORD。

Answer 2

我遇到了这个：http://arma.sourceforge.net/docs.html#config_hpp。

解决方法是在犰狳的头文件上方添加一行来设置64位整数选项：

#define ARMA_64BIT_WORD

我把它放在 [R root]/lib/R/library/RcppArmadillo/include/RcppArmadilloConfig.h 中

我曾尝试为 cxxfunction 使用“includes”参数，但我认为这个定义语句必须在 include 语句之上，

#include RcppArmadillo.h

在 cpp 代码中。我不知道内联包的cxxfunction是否允许这样做。

修改 RcppArmadilloConfig.h 后，我现在可以声明一个大于 2^16 的矩阵。

> m <- 1000
> n <- 65536+1000
> 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:
10218.8
NULL
>