Matrix Multiplication Using NumericMatrix and NumericVector in Rcpp

Question

I am wondering is there a way of calculating matrix multiplication using NumericMatrix and NumericVector class. I am wondering if there is any simple way to help me avoid the following loop to conduct this calculation. I just want to calculate X%*%beta.

// assume X and beta are initialized and X is of dimension (nsites, p), 
// beta is a NumericVector with p elements. 
for(int j = 0; j < nsites; j++)
 {
    temp = 0;

    for(int l = 0; l < p; l++) temp = temp + X(j,l) * beta[l];

}

Thank you very much in advance!

I would look into [RcppArmadillo](http://dirk.eddelbuettel.com/code/rcpp.armadillo.html) or RcppEigen. — Dirk Eddelbuettel, Feb 11 '15 at 22:31
I see, just to confirm, Rcpp sugar does not have %*% like R, right? Thank you a lot for your help! — Crystal, Feb 11 '15 at 22:32

score 5 · Accepted Answer · answered Feb 12 '15 at 15:55

Building off of Dirk's comment, here are a few cases that demonstrate the Armadillo library's matrix multiplication via the overloaded * operator:

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

// [[Rcpp::export(".mm")]]
arma::mat mm_mult(const arma::mat& lhs,
                  const arma::mat& rhs)
{
  return lhs * rhs;
}

// [[Rcpp::export(".vm")]]
arma::mat vm_mult(const arma::vec& lhs,
                  const arma::mat& rhs)
{
  return lhs.t() * rhs;
}

// [[Rcpp::export(".mv")]]
arma::mat mv_mult(const arma::mat& lhs,
                  const arma::vec& rhs)
{
  return lhs * rhs;
}

// [[Rcpp::export(".vv")]]
arma::mat vv_mult(const arma::vec& lhs,
                  const arma::vec& rhs)
{
  return lhs.t() * rhs;
}

You could then define an R function to dispatch the appropriate C++ function:

`%a*%` <- function(x,y) {

  if (is.matrix(x) && is.matrix(y)) {
    return(.mm(x,y))
  } else if (!is.matrix(x) && is.matrix(y)) {
    return(.vm(x,y))
  } else if (is.matrix(x) && !is.matrix(y)) {
    return(.mv(x,y))
  } else {
    return(.vv(x,y))
  }

}
##
mx <- matrix(1,nrow=3,ncol=3)
vx <- rep(1,3)
my <- matrix(.5,nrow=3,ncol=3)
vy <- rep(.5,3)

And comparing to R's %*% function:

R>  mx %a*% my
     [,1] [,2] [,3]
[1,]  1.5  1.5  1.5
[2,]  1.5  1.5  1.5
[3,]  1.5  1.5  1.5

R>  mx %*% my
     [,1] [,2] [,3]
[1,]  1.5  1.5  1.5
[2,]  1.5  1.5  1.5
[3,]  1.5  1.5  1.5
##
R>  vx %a*% my
     [,1] [,2] [,3]
[1,]  1.5  1.5  1.5

R>  vx %*% my
     [,1] [,2] [,3]
[1,]  1.5  1.5  1.5
##
R>  mx %a*% vy
     [,1]
[1,]  1.5
[2,]  1.5
[3,]  1.5

R>  mx %*% vy
     [,1]
[1,]  1.5
[2,]  1.5
[3,]  1.5
##
R>  vx %a*% vy
     [,1]
[1,]  1.5

R>  vx %*% vy
     [,1]
[1,]  1.5

Thank you so much! This demonstration is very helpful and clear to a beginner like me! I highly appreciate your help! — Crystal, Feb 13 '15 at 03:24
Because there was no intention of modifying the arguments. Passing by `const&` is a common idiom, see any number of discussions such as [here](http://stackoverflow.com/questions/5060137/passing-as-const-and-by-reference-worth-it) and [here](http://stackoverflow.com/questions/136880/sell-me-on-const-correctness). — nrussell, Oct 06 '16 at 21:24

Matrix Multiplication Using NumericMatrix and NumericVector in Rcpp

1 Answers1

Linked

Related