犰狳:从稀疏矩阵中获取稀疏行向量的非零位置

Question

我正在通过RcppArmadillo使用犰狳和R。

我有一个稀疏矩阵和一个行号作为输入。我想搜索矩阵的相应行并返回非零的位置。

到目前为止，我的函数看起来像

// [[Rcpp::export]]
arma::uvec findAdjacentStates(arma::sp_mat adjacency, int row) {
  arma::uvec out(adjacency.n_cols);
  arma::SpSubview<double>::const_iterator start = adjacency.row(row).begin();
  arma::SpSubview<double>::const_iterator end = adjacency.row(row).end();
  for(arma::SpSubview<double>::const_iterator it = start; it != end; ++it)
  {
    Rcout << "location: " << it.row() << "," << it.col() << "  ";
    Rcout << "value: " << (*it) << std::endl;
  }
  return out;
}

这是基于之前的 SO 答案。

该函数使 R 崩溃。

require(Matrix)
x = rsparsematrix(10, 10, .2)
x = x > 1
x = as(x, "dgCMatrix")
findAdjacentStates(x, 1)

我不清楚的一件事是如何专门迭代行向量;之前的 SO 答案是迭代稀疏矩阵。

更新：根据 valgrind 的说法，问题出在 operator++ （SpSubview_iterators_meat.hpp：319） 中，所以这似乎不是迭代稀疏行向量的正确方法

Answer 1

迭代稀疏矩阵行的方法是使用 sp_mat：：row_iterator。不幸的是，无法提前知道输出向量的大小，并且uvec对象不像常规向量那样具有push_back。这是我的建议：

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
using namespace arma;
// [[Rcpp::export]]
IntegerVector findAdjacentStates(sp_mat adjacency, int row) {
    IntegerVector out;
    sp_mat::const_row_iterator start = adjacency.begin_row(row);
    sp_mat::const_row_iterator end = adjacency.end_row(row);
    for ( sp_mat::const_row_iterator i = start; i != end; ++i )
    {
        out.push_back(i.col());
    }
    return out;
}

我们可以很容易地测试出来：

# We need Rcpp and Matrix:
library(Rcpp)
library(Matrix)
# This is the C++ script I wrote:
sourceCpp('SOans.cpp')
# Make example data (setting seed for reproducibility):
set.seed(123)
x = rsparsematrix(10, 10, .2)
# And test the function:
findAdjacentStates(x, 1)
[1] 4
x
10 x 10 sparse Matrix of class "dgCMatrix"
 [1,]  .    .  0.84 .  0.40  0.7 .  .     .    -0.56
 [2,]  .    .  .    . -0.47  .   .  .     .     .   
 [3,]  .    .  .    .  .     .   .  .    -2.00  .   
 [4,]  0.15 .  .    .  .     .   .  .     .    -0.73
 [5,]  1.80 .  .    .  .     .   .  .     .     .   
 [6,]  .    .  .    .  .     .   .  .     0.11  .   
 [7,]  .    . -1.10 .  .     .   . -1.70 -1.10  .   
 [8,]  .    .  .    .  .     .   .  1.30  .    -0.22
 [9,] -0.63 .  1.20 .  .     .   .  0.36  .     .   
[10,]  .    .  .    .  0.50 -1.0 .  .     .     .   

所以，我们可以看到这是有效的;第 1 行（C++ 术语;第 2 行 R 术语）只有一个非零元素，在第 4 列（C++ 术语中;第 5 列在 R 术语中）。如果要将输出返回到 R，这应该有效。如果您想在另一个C++函数中使用输出，根据您正在执行的操作，您可能更喜欢使用uvec而不是IntegerVector，但您可以将IntegerVector转换为uvec（可能不是最有效的解决方案，但我现在想到的最好的解决方案）：

// [[Rcpp::export]]
uvec findAdjacentStates(sp_mat adjacency, int row) {
    IntegerVector tmp;
    sp_mat::const_row_iterator start = adjacency.begin_row(row);
    sp_mat::const_row_iterator end = adjacency.end_row(row);
    for ( sp_mat::const_row_iterator i = start; i != end; ++i )
    {
        tmp.push_back(i.col());
    }
    uvec out = as<uvec>(tmp);
    return out;
}