代码之家 › 专栏 › 技术社区 › JCWong

Rcpp特征稀疏矩阵cbind

rcpp eigen matrix r

JCWong · 技术社区 · 8 年前

link ,但它适用于稠密矩阵上的cbind

#include <RcppEigen.h>
// [[Rcpp::depends(RcppEigen)]]

using namespace Rcpp;
using namespace Eigen;

// [[Rcpp::export]]
Eigen::SparseMatrix<double> RcppMatrixCbind(Eigen::MappedSparseMatrix<double>& matrix1,
                                            Eigen::MappedSparseMatrix<double>& matrix2) {

  SparseMatrix<double> out(matrix1.rows(), matrix1.cols() + matrix2.cols());
  std::vector<Triplet<double> > tripletList;
  tripletList.reserve(matrix1.nonZeros() + matrix2.nonZeros());
  for (int k = 0; k < matrix1.outerSize(); ++k)
  {
    for (MappedSparseMatrix<double>::InnerIterator it(matrix1, k); it; ++it)
    {
      tripletList.push_back(Triplet<double>(it.row(), it.col(), it.value()));
    }
  }
  for (int k = 0; k < matrix2.outerSize(); ++k)
  {
    for (MappedSparseMatrix<double>::InnerIterator it(matrix2, k); it; ++it)
    {
      tripletList.push_back(Triplet<double>(it.row(), it.col() + matrix1.cols(), it.value()));
    }
  }
  out.setFromTriplets(tripletList.begin(), tripletList.end());
  return out;
}

/*** R
require(Matrix)
x = rsparsematrix(10000, 500, .1)
system.time(foo <- cbind(x,x))
system.time(bar <- RcppMatrixCbind(x, x))
  */

1 回复 | 直到 8 年前

Dirk is no longer here 8 年前

这实际上更像是一个特征问题:如何扩展稀疏矩阵?

您在解决方案中所做的是逐元素执行所有操作,很可能是可以打败它。这就是我们看到的 Matrix 解决方案,转到高效代码,可能来自CHOLMD。

我快速查看了Eigen文档。它警告说:

关于读取访问,稀疏矩阵公开了与密集矩阵相同的API来访问子矩阵,如块、列和行。有关详细介绍,请参阅块操作。然而,出于性能原因,写入子稀疏矩阵受到更大的限制,目前只有列主稀疏矩阵(行主稀疏矩阵)的连续列集(行)是可写的。此外,这些信息必须在编译时已知,而不包括以下方法: block(...) 和 corner*(...)

cbind() 在blocky够了。一个非常简单的解决方案是

// [[Rcpp::export]]
Eigen::SparseMatrix<double>
RcppMatrixCb2(Eigen::MappedSparseMatrix<double>& matrix1,
              Eigen::MappedSparseMatrix<double>& matrix2) {

  SparseMatrix<double> out(matrix1.rows(), matrix1.cols() + matrix2.cols());

  out.leftCols(matrix1.cols()) = matrix1; 
  out.rightCols(matrix2.cols()) = matrix2; 
  return out;
}

> benchmark(Matrix=cbind(x,x), 
+           prevSol=RcppMatrixCbind(x,x), 
+           newSol=RcppMatrixCb2(x,x),
+           order="relative")[,1:4]
     test replications elapsed relative
3  newSol          100   0.585    1.000
1  Matrix          100   0.686    1.173
2 prevSol          100   4.603    7.868
>

我的完整文件如下。它在两个函数中都缺乏测试:我们需要确保矩阵1和矩阵2具有相同的行。

// cf https://stackoverflow.com/questions/45875668/rcpp-eigen-sparse-matrix-cbind

#include <RcppEigen.h>
// [[Rcpp::depends(RcppEigen)]]

using namespace Rcpp;
using namespace Eigen;

// [[Rcpp::export]]
Eigen::SparseMatrix<double> RcppMatrixCbind(Eigen::MappedSparseMatrix<double>& matrix1,
                                            Eigen::MappedSparseMatrix<double>& matrix2) {

  SparseMatrix<double> out(matrix1.rows(), matrix1.cols() + matrix2.cols());
  std::vector<Triplet<double> > tripletList;
  tripletList.reserve(matrix1.nonZeros() + matrix2.nonZeros());
  for (int k = 0; k < matrix1.outerSize(); ++k) {
    for (MappedSparseMatrix<double>::InnerIterator it(matrix1, k); it; ++it) {
      tripletList.push_back(Triplet<double>(it.row(), it.col(), it.value()));
    }
  }
  for (int k = 0; k < matrix2.outerSize(); ++k) {
    for (MappedSparseMatrix<double>::InnerIterator it(matrix2, k); it; ++it) {
      tripletList.push_back(Triplet<double>(it.row(), it.col() + matrix1.cols(), it.value()));
    }
  }
  out.setFromTriplets(tripletList.begin(), tripletList.end());
  return out;
}

// [[Rcpp::export]]
Eigen::SparseMatrix<double> RcppMatrixCb2(Eigen::MappedSparseMatrix<double>& matrix1,
                                          Eigen::MappedSparseMatrix<double>& matrix2) {

  SparseMatrix<double> out(matrix1.rows(), matrix1.cols() + matrix2.cols());

  out.leftCols(matrix1.cols()) = matrix1; 
  out.rightCols(matrix2.cols()) = matrix2; 
  return out;
}


/*** R
require(Matrix)
set.seed(42)
x = rsparsematrix(10000, 500, .1)
library(rbenchmark)
benchmark(Matrix=cbind(x,x), 
          prevSol=RcppMatrixCbind(x,x), 
          newSol=RcppMatrixCb2(x,x),
          order="relative")[,1:4]
*/