Матричное умножение с использованием NumericMatrix и NumericVector в Rcpp

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]; 

} 

спасибо, что заблаговременно!

Answer 1

здание от комментария Дирка, вот несколько случаев, которые демонстрируют умножение матриц библиотеки 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> 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