列表示 3d 矩阵（立方体）Rcpp答案

【问题标题】：Column means 3d matrix (cube) Rcpp列表示 3d 矩阵（立方体）Rcpp
【发布时间】：2015-06-21 22:07:12
【问题描述】：

我有一个程序，我需要在 Rcpp 中重复计算立方体X(nRow, nCol, nSlice) 的每个切片的列均值，结果均值形成矩阵M(nCol, nSlice)。以下代码产生了错误：

#include <RcppArmadillo.h>

// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp; 
using namespace arma;

// [[Rcpp::export]]

mat cubeMeans(arma::cube X){
   int nSlice = X.n_slices;
   int nCol = X.n_cols;
   int nRow = X.n_rows;
   arma::vec Vtmp(nCol);
   arma::mat Mtmp(nRow, nCol);
   arma::mat Means(nCol, nSlice);
   for (int i = 0; i < nSlice; i++){
      Mtmp = X.slice(i);
      for(int j = 0; j < nCol; j++){
         Vtmp(j) = sum(Mtmp.col(j))/nRow; 
      }
      Means.col(i) = Vtmp;
   }
  return(wrap(Means));
}

'/Rcpp/internal/Exporter.h:31:31: 错误：没有匹配函数调用'arma::Cube::Cube(SEXPREC*&)'

我想不通。当函数的输入是矩阵（并返回向量）时，我没有收到错误。但是，我将上述函数作为我的主程序的一部分，即

#include <RcppArmadillo.h>

// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
using namespace arma;

mat cubeMeans(arma::cube X){
  int nSlice = X.n_slices;
  ...
  return(Means);
}

// [[Rcpp::export]]

main part of program

程序编译成功，但是速度很慢（几乎和使用colMeans的R版本程序一样慢）。有没有更好的方法来计算多维数据集上的列均值，为什么会出现编译错误？

我将不胜感激。

问候，

【问题讨论】：

几乎和 colMeans 一样慢？？？ ……那就这样吧。 colMeans 在 C 中进行了优化。您没有提供用于运行该程序的代码可能包含回答问题的线索。顺便说一句，R 中没有“3d 矩阵”。数组将是另一回事。术语通常很重要。

标签： r matrix rcpp armadillo

【解决方案1】：

我在尝试使用 arma::cube 作为 Rcpp 函数参数时也收到此错误。 ~~基于编译器错误，我相信这是因为当前没有定义Rcpp::wrap<arma::cube>（需要它来处理您将传递给函数的 R 对象）。~~† 在阅读了一些相关的之后在线示例，看起来典型的解决方法是将您的 R array 作为NumericVector 读取，并且由于它保留其dims 属性，因此使用这些来设置您的arma::cube 尺寸。尽管需要一两个额外的步骤来解决 ~~缺少的 wrap 专业化~~†，但我组合的 Armadillo 版本似乎比我的 R 解决方案快很多：

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

// [[Rcpp::export]]
arma::mat cube_means(Rcpp::NumericVector vx) {

  Rcpp::IntegerVector x_dims = vx.attr("dim");
  arma::cube x(vx.begin(), x_dims[0], x_dims[1], x_dims[2], false);

  arma::mat result(x.n_cols, x.n_slices);
  for (unsigned int i = 0; i < x.n_slices; i++) {
    result.col(i) = arma::conv_to<arma::colvec>::from(arma::mean(x.slice(i)));  
  }

  return result;
}

/*** R

rcube_means <- function(x) t(apply(x, 2, colMeans))

xl <- array(1:10e4, c(100, 100 ,10))
all.equal(rcube_means(xl), cube_means(xl))
#[1] TRUE

R> microbenchmark::microbenchmark(
    "R Cube Means" = rcube_means(xl),
    "Arma Cube Means" = cube_means(xl),
    times = 200L)
Unit: microseconds
            expr      min       lq      mean   median       uq       max neval
    R Cube Means 6856.691 8204.334 9843.7455 8886.408 9859.385 97857.999   200
 Arma Cube Means  325.499  380.540  643.7565  416.863  459.800  3068.367   200

*/

我利用arma::mats 的arma::mean 函数重载将默认计算列均值这一事实（arma::mean(x.slice(i), 1) 将为您提供该切片的行均值）。

编辑： † 再想一想，我不确定这是否与 Rcpp::wrap 有关 - 但问题似乎与缺少 Exporter<> 专业化有关arma::cube - Rcpp 的 Exporter.h 第 31 行：

template <typename T>
class Exporter{
public:
  Exporter( SEXP x ) : t(x){}
  inline T get(){ return t ; }

private:
  T t ;
} ;

不管怎样，NumericVector / 我使用的设置尺寸方法目前似乎是功能性解决方案。

根据您在问题中描述的输出维度，我假设您希望结果矩阵的每一列都是相应数组切片的列均值向量（列 1 = 切片 1 的列均值等... )，即

R> x <- array(1:27, c(3, 3, 3))
R> rcube_means(x)
     [,1] [,2] [,3]
[1,]    2   11   20
[2,]    5   14   23
[3,]    8   17   26
R> cube_means(x)
     [,1] [,2] [,3]
[1,]    2   11   20
[2,]    5   14   23
[3,]    8   17   26

但如果需要，您可以轻松更改它。

【讨论】：