我在想,是否有一種使用NumericMatrix和NumericVector類來計算矩陣乘法的方法。我想知道是否有任何簡單的方法 來幫助我避免以下循環來執行此計算。我只是想計算X%*%beta。矩陣乘法在Rcpp中使用NumericMatrix和NumericVector
非常感謝您提前!
我在想,是否有一種使用NumericMatrix和NumericVector類來計算矩陣乘法的方法。我想知道是否有任何簡單的方法 來幫助我避免以下循環來執行此計算。我只是想計算X%*%beta。矩陣乘法在Rcpp中使用NumericMatrix和NumericVector
非常感謝您提前!
大廈關閉德克的評論,在這裏是通過重載*
運營商展示了犰狳庫的矩陣乘法少數病例:
#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
非常感謝你!這個演示對我這樣的初學者非常有幫助和清晰!我非常感謝你的幫助! – Crystal 2015-02-13 03:24:22
爲什麼使用'const'? – Jason 2016-10-06 19:59:50
因爲沒有修改參數的意圖。通過'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
我會去了解一下[RcppArmadillo(http://dirk.eddelbuettel.com/code/rcpp.armadillo.html)或RcppEigen。 – 2015-02-11 22:31:23
我看到了,只是爲了確認,Rcpp糖沒有像R那樣的%*%,對吧?非常感謝您的幫助! – Crystal 2015-02-11 22:32:55