在Rcpp中使用NumericMatrix和NumericVector进行矩阵乘法

Question

我想知道是否有一种使用NumericMatrix和NumericVector类计算矩阵乘法的方法。 是否有简单的方法可以帮助我避免以下循环以进行此计算。 我只想计算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];

}

非常感谢您的提前帮助！以下是需要翻译的内容：

谢谢！

- Crystal

---

我会研究RcppArmadillo或RcppEigen。 - Dirk Eddelbuettel

我明白了，只是确认一下，Rcpp sugar 没有像 R 一样的 %*% 运算符，对吗？非常感谢您的帮助！ - Crystal

Answer 1

在 Dirk 的评论基础上，以下是几个案例，演示了 Armadillo 库通过重载的 * 运算符进行矩阵乘法：

#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;
}

您可以定义一个R函数来分派相应的C++函数：

`%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)

而与R的%*%函数相比：

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

- nrussell

---

非常感谢！这个演示对像我这样的初学者非常有帮助和清晰！我非常感激你的帮助！ - Crystal

为什么要使用 const？ - Jason

因为没有修改参数的意图，所以通过const&传递是一种常见的习惯用法，请参考许多讨论，例如这里和这里。 - nrussell