Search code examples
rrcpp

Insert submatrix with Rcpp


I am trying to achieve the following R example with Rcpp:

X <- matrix(0, 5, 10)
X[1:4, 4] <- rexp(4)

What I have tried so far is:

#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericMatrix foo1() {
  NumericMatrix X(5, 10);
  NumericMatrix y(4, 1);
  y(_, 0) = rexp(4, 1);
  X(Range(0,3),Range(3,3)) = y;
  return X; 
}

But I keep getting a compilation error, saying that no match for 'operator='. Any hints to what I am doing wrong would be greatly appreciated.


Solution

  • Matrix operations with Rcpp are a bit lacking for better or worse. Any in-depth matrix work should be done with either RcppArmadillo or RcppEigen.

    Sample implementation:

    #include <RcppArmadillo.h>
    // [[Rcpp::depends(RcppArmadillo)]]
    
    // [[Rcpp::export]]
    arma::mat matrix_fill_single_col() {
    
      // Setup X matrix
      arma::mat X = arma::zeros<arma::mat>(5, 10);
    
      // Generate random values from exponential and save into a vector.
      arma::vec y = Rcpp::as<arma::vec>(Rcpp::rexp(4, 1));
    
      // Fill the fourth column in X (Recall: C++ indexes start at 0 not 1)
      X.submat(0, 3, 3, 3) = y;
      // Or...
      // X.col(3) = y;
    
      return X; 
    }
    

    Test

    matrix_fill_single_col()
    #      [,1] [,2] [,3]      [,4] [,5] [,6] [,7] [,8] [,9] [,10]
    # [1,]    0    0    0 0.2685970    0    0    0    0    0     0
    # [2,]    0    0    0 1.6018346    0    0    0    0    0     0
    # [3,]    0    0    0 0.6467853    0    0    0    0    0     0
    # [4,]    0    0    0 0.6655340    0    0    0    0    0     0
    # [5,]    0    0    0 0.0000000    0    0    0    0    0     0