Sparse double matrix times vector of vars

Developers
1

This works with RcppEigen 3.3.x. We can talk about generalizing it for rstanarm at the meeting.

X <- rbind(
  c(0,	3,	0,	0,	0),
  c(22,	0,	0,	0,	17),
  c(7,	5,	0,	1,	0),
  c(0,	0,	0,	0,	0),
  c(0,	0,	14,	0,	8)
)
x <- c(X)
x <- x[x != 0]

InnerIndices <- c(1L,	2L,	0L,	2L,	4L,	2L,	1L,	4L)
OuterStarts <- c(0L,	2L,	4L,	5L,	6L,	8L)

beta <- rnorm(nrow(X))
Rcpp::sourceCpp('sparse.cpp')
args(smv)
cbind(R = X %*% beta, Eigen = 
        smv(rows = nrow(X), cols = ncol(X), Values = x, 
            InnerIndices = InnerIndices, OuterStarts = OuterStarts, 
            InnerNNZs = length(x), beta = beta) )

where sparse.cpp is a file that contains

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

Eigen::VectorXd
smv(int rows, int cols,
    std::vector<double> Values, std::vector<int> InnerIndices, std::vector<int> OuterStarts,
    int InnerNNZs, Eigen::VectorXd beta) {
  
  const Eigen::Index r = rows;
  const Eigen::Index c = cols;
  const Eigen::Index nnz = InnerNNZs;

  double* Values_ptr = &Values[0];
  typedef Eigen::SparseMatrix<double>::StorageIndex StorageIndex;
  StorageIndex* OS_ptr = &OuterStarts[0];
  StorageIndex* II_ptr = &InnerIndices[0];
  
  Eigen::Map<const Eigen::SparseMatrix<double> > 
  sm1(r, c, nnz, OS_ptr, II_ptr, Values_ptr);
  return sm1 * beta;
}
2

Here is an improved version of the C++ code

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

Eigen::VectorXd
smv(int rows, int cols,
    const std::vector<double>& Values,
    const std::vector<int>& InnerIndices, 
    const std::vector<int>& OuterStarts, 
    const Eigen::VectorXd& beta) {
  
  Eigen::Map<const Eigen::SparseMatrix<double,Eigen::RowMajor> > 
  sm(rows, cols, Values.size(), &OuterStarts[0], &InnerIndices[0], &Values[0]);

  // To do derivatives, iterate over the non-zeros with this pattern from
  // http://eigen.tuxfamily.org/dox/group__TutorialSparse.html
  for (int k = 0; k < sm.outerSize(); k++) {
    for (Eigen::Map<const Eigen::SparseMatrix<double, Eigen::RowMajor> >::InnerIterator 
           it(sm,k); it; ++it) {
      std::cout << "Derivative of output " << k << 
        " w.r.t. input " << it.col() << " = " << it.value() << std::endl;
    }
  }

  return sm * beta;
}

and you can call it from R like

X <- rbind(
  c(0,	3,	0,	0,	0),
  c(22,	0,	0,	0,	17),
  c(7,	5,	0,	1,	0),
  c(0,	0,	0,	0,	0),
  c(0,	0,	14,	0,	8)
)
beta <- rnorm(nrow(X))
parts <- rstan::extract_sparse_parts(X)

Rcpp::sourceCpp('sparse.cpp')
args(smv)
cbind(R = c(X %*% beta), Eigen = 
        smv(rows = nrow(X), cols = ncol(X), Values = parts$w, 
            # Remember that Eigen starts counting at zero
            InnerIndices = parts$v - 1L, OuterStarts = parts$u - 1L, 
            beta = beta) )
1 Like