将 Matrix::sparseMatrix 传递给 Rcpp
Posted
技术标签:
【中文标题】将 Matrix::sparseMatrix 传递给 Rcpp【英文标题】：passing Matrix::sparseMatrix to Rcpp 【发布时间】：2020-03-24 20:40:05 【问题描述】：所以我对将稀疏矩阵从 R 传递到 c++ 的推荐方法感到非常困惑。我的印象是 sp_mat 是正确的参数类型，如下面的代码所示
testCode = '
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
void testFun(arma::sp_mat F)
Rcpp::Rcout << "F has " << F.n_rows << " rows" << std::endl;
'
Rcpp::sourceCpp(code = testCode)
n = 70000
M = Matrix::sparseMatrix(i=c(n), j=c(n), x=c(1))
testFun(M)
但是，运行此代码会产生以下错误：
error: SpMat::init(): requested size is too large
Error in testFun(M) : SpMat::init(): requested size is too large
Calls: testFun -> .Call
Execution halted
我在https://gallery.rcpp.org/articles/armadillo-sparse-matrix/ 中看到了这个例子，但我不确定它是否是说每次我们将稀疏矩阵传递给 c++ 时，我们都应该使用那里提供的函数？感谢您的澄清！
【问题讨论】：
您遇到了运行时错误。 N=70000 不通过。尝试使用较小的 N。查看现有示例以了解如何使用它 - 请注意，例如，您提供的链接具有 索引向量 作为i 和 j 的参数。
我的意思是我需要我的矩阵有 N 这么大（甚至实际上更大）。让我用向量来改变这个例子，但这真的不是问题。我实际需要传递的矩阵是自动创建的“dgCMatrix”实例。
一个 dense 矩阵将有行 * 列 * 字节。 sparse 矩阵本质上具有 3 * N * 8，因为它将索引向量与单元格值一起保存。因此，标称大小可以为 70k x 70k，并且如果您实际提供单元格索引作为稀疏格式需要，仍然适合它。
是的，没错！这正是我在这里使用稀疏矩阵的原因！此示例中的矩阵 M 的尺寸为 70k x 70k，但只有一个非零值 - 这就是我对错误感到惊讶的原因。看起来我不理解某些东西，因为在我看来，我的示例（现在我将索引更改为向量 - 没有帮助）与 Rcpp 画廊中发布的示例相同，除了 nominal i> 矩阵的维度。
好的，我想我在犰狳源代码中找到了答案。原来这里也有描述：***.com/questions/40592054/…
【参考方案1】：
好的，所以我想我找到了答案。基本上，如果元素的总数大于存储元素数量的变量的大小，则犰狳会抛出此错误，如下所示： https://gitlab.com/conradsnicta/armadillo-code/-/blob/9.900.x/include/armadillo_bits/SpMat_meat.hpp
之前有人意识到这一点，并在此处提供了解决方案： Large Matrices in RcppArmadillo via the ARMA_64BIT_WORD define
【讨论】：
【参考方案2】：如果您重新访问基本的example in the Rcpp Gallery 并设置一个或两个稀疏矩阵对象，很明显j 的高值会导致p 槽中的完全扩展（检查从 @ 返回的对象987654324@).
所以这是一个更简单的例子，它的i 值（仍然很高）但j 的值很低。我想你应该可以从这里拿走它：
代码
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp ;
// [[Rcpp::export]]
void convertSparse(S4 mat)
// obtain dim, i, p. x from S4 object
IntegerVector dims = mat.slot("Dim");
arma::urowvec i = Rcpp::as<arma::urowvec>(mat.slot("i"));
arma::urowvec p = Rcpp::as<arma::urowvec>(mat.slot("p"));
arma::vec x = Rcpp::as<arma::vec>(mat.slot("x"));
int nrow = dims[0], ncol = dims[1];
// use Armadillo sparse matrix constructor
arma::sp_mat res(i, p, x, nrow, ncol);
Rcout << "SpMat res:\n" << res << std::endl;
// [[Rcpp::export]]
void convertSparse2(S4 mat) // slight improvement with two non-nested loops
IntegerVector dims = mat.slot("Dim");
arma::urowvec i = Rcpp::as<arma::urowvec>(mat.slot("i"));
arma::urowvec p = Rcpp::as<arma::urowvec>(mat.slot("p"));
arma::vec x = Rcpp::as<arma::vec>(mat.slot("x"));
int nrow = dims[0], ncol = dims[1];
arma::sp_mat res(nrow, ncol);
// create space for values, and copy
arma::access::rw(res.values) = arma::memory::acquire_chunked<double>(x.size() + 1);
arma::arrayops::copy(arma::access::rwp(res.values), x.begin(), x.size() + 1);
// create space for row_indices, and copy
arma::access::rw(res.row_indices) = arma::memory::acquire_chunked<arma::uword>(i.size() + 1);
arma::arrayops::copy(arma::access::rwp(res.row_indices), i.begin(), i.size() + 1);
// create space for col_ptrs, and copy
arma::access::rw(res.col_ptrs) = arma::memory::acquire<arma::uword>(p.size() + 2);
arma::arrayops::copy(arma::access::rwp(res.col_ptrs), p.begin(), p.size() + 1);
// important: set the sentinel as well
arma::access::rwp(res.col_ptrs)[p.size()+1] = std::numeric_limits<arma::uword>::max();
// set the number of non-zero elements
arma::access::rw(res.n_nonzero) = x.size();
Rcout << "SpMat res:\n" << res << std::endl;
/*** R
suppressMessages(
library(methods)
library(Matrix)
)
i <- c(1,3:6)
j <- c(2,9,6:8)
x <- 5 * (1:5)
A <- sparseMatrix(i, j, x = x)
print(A)
convertSparse(A)
i <- 56789
j <- 87
x <- 42
B <- sparseMatrix(i, j, x=x)
#print(B)
convertSparse(B)
convertSparse2(B)
*/
输出
R> Rcpp::sourceCpp("~/git/***/60838958/answer.cpp")
R> suppressMessages(library(methods); library(Matrix))
R> i <- c(1,3:6)
R> j <- c(2,9,6:8)
R> x <- 5 * (1:5)
R> A <- sparseMatrix(i, j, x = x)
R> print(A)
6 x 9 sparse Matrix of class "dgCMatrix"
[1,] . 5 . . . . . . .
[2,] . . . . . . . . .
[3,] . . . . . . . . 10
[4,] . . . . . 15 . . .
[5,] . . . . . . 20 . .
[6,] . . . . . . . 25 .
R> convertSparse(A)
SpMat res:
[matrix size: 6x9; n_nonzero: 5; density: 9.26%]
(0, 1) 5.0000
(3, 5) 15.0000
(4, 6) 20.0000
(5, 7) 25.0000
(2, 8) 10.0000
R> i <- 56789
R> j <- 87
R> x <- 42
R> B <- sparseMatrix(i, j, x=x)
R> #print(B)
R> convertSparse(B)
SpMat res:
[matrix size: 56789x87; n_nonzero: 1; density: 2.02e-05%]
(56788, 86) 42.0000
R> convertSparse2(B)
SpMat res:
[matrix size: 56789x87; n_nonzero: 1; density: 2.02e-05%]
(56788, 86) 42.0000
R>
编辑确实，很好的提醒。如果我们添加
#define ARMA_64BIT_WORD 1
在包含 RcppArmadillo.h 标头之前，然后使用 i 和 j large 一切正常。下面输出的尾部。
更新的输出（只是尾端）
R> i <- 56789
R> j <- 87654
R> x <- 42
R> B <- sparseMatrix(i, j, x=x)
R> #print(B)
R> convertSparse(B)
SpMat res:
[matrix size: 56789x87654; n_nonzero: 1; density: 2.01e-08%]
(56788, 87653) 42.0000
R> convertSparse2(B)
SpMat res:
[matrix size: 56789x87654; n_nonzero: 1; density: 2.01e-08%]
(56788, 87653) 42.0000
R>
【讨论】：
以上是关于将 Matrix::sparseMatrix 传递给 Rcpp的主要内容，如果未能解决你的问题，请参考以下文章
为啥将 ClientID 传递给 javascript 函数会传递整个控件？
将对象从 Activity 传递到 Fragments 是通过引用传递
如何将列表的所有元素作为函数的参数传递而不将列表作为参数传递