2015-02-11 96 views
3

我在想,是否有一种使用NumericMatrix和NumericVector类来计算矩阵乘法的方法。我想知道是否有任何简单的方法 来帮助我避免以下循环来执行此计算。我只是想计算X%*%beta。矩阵乘法在Rcpp中使用NumericMatrix和NumericVector

​​

非常感谢您提前!

+0

我会去了解一下[RcppArmadillo(http://dirk.eddelbuettel.com/code/rcpp.armadillo.html)或RcppEigen。 – 2015-02-11 22:31:23

+0

我看到了,只是为了确认,Rcpp糖没有像R那样的%*%,对吧?非常感谢您的帮助! – Crystal 2015-02-11 22:32:55

回答

3

大厦关闭德克的评论,在这里是通过重载*运营商展示了犰狳库的矩阵乘法少数病例:

#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 
+1

非常感谢你!这个演示对我这样的初学者非常有帮助和清晰!我非常感谢你的帮助! – Crystal 2015-02-13 03:24:22

+0

为什么使用'const'? – Jason 2016-10-06 19:59:50

+0

因为没有修改参数的意图。通过'const&'传递是一个常见的习惯用法,请参阅任何数量的讨论,如[here](http://stackoverflow.com/questions/5060137/passing-as-const-and-by-reference-worth-it)和[这里](http://stackoverflow.com/questions/136880/sell-me-on-const-correctness)。 – nrussell 2016-10-06 21:24:22