使用 RcppArmadillo 在矩阵的列上应用函数有效，但在应用于行时返回错误

Posted 2023-02-16

【中文标题】使用 RcppArmadillo 在矩阵的列上应用函数有效，但在应用于行时返回错误

【英文标题】：Applying a function over columns of matrix using RcppArmadillo works but returns error when applied over rows

【发布时间】：2019-02-24 13:28:32

【问题描述】：

我在Rcpp 中编写了一个函数qSelectMbycol，它在O(n) 时间内返回每列的kth 最大元素。此功能工作正常。如果我尝试做同样的事情但处理行而不是列（函数qSelectMbyrow），它会返回错误"error: Mat::init(): requested size is not compatible with column vector layout"。有人想我做错了什么吗？我将此文件保存为“qselect.cpp”：

// [[Rcpp::depends(RcppArmadillo)]]
#define RCPP_ARMADILLO_RETURN_COLVEC_AS_VECTOR
#include <RcppArmadillo.h>
using namespace arma;

// [[Rcpp::export]]
arma::vec qSelectMbycol(arma::mat& M, const int k) 

  // ARGUMENTS
  // M: matrix for which we want to find the k-th largest elements of each column
  // k: k-th statistic to look up

  arma::mat Y(M.memptr(), M.n_rows, M.n_cols);
  // we apply over columns
  int c = M.n_cols;
  arma::vec out(c);
  int i;
  for (i = 0; i < c; i++) 
      arma::vec y = Y.col(i);
      std::nth_element(y.begin(), y.begin() + k - 1, y.end());
      out(i) = y(k-1); // the k-th largest value of each column
  

  return out;



// [[Rcpp::export]]
arma::vec qSelectMbyrow(arma::mat& M, const int k) 

  // ARGUMENTS
  // M: matrix for which we want to find the k-th largest elements of each row
  // k: k-th statistic to look up

  arma::mat Y(M.memptr(), M.n_rows, M.n_cols);
  // we apply over rows
  int r = M.n_rows;
  arma::vec out(r);
  int i;
  for (i = 0; i < r; i++) 
    arma::vec y = Y.row(i); // this line throws the error "error: Mat::init(): requested size is not compatible with column vector layout"
    std::nth_element(y.begin(), y.begin() + k - 1, y.end());
    out(i) = y(k-1); // should give k-th largest value of each row
  

  return out;

例子：

n=500
p=100
set.seed(1)
M=matrix(rnorm(n, mean = 100, sd = 1),n,1)
library(Rcpp)
library(RcppArmadillo)
Rcpp::sourceCpp('qselect.cpp')
qSelectMbycol(M,5) # works OK
qSelectMbyrow(M,5) # throws error "error: Mat::init(): requested size is not compatible with column vector layout"

我也试过插入

  typedef std::vector<double> stdvec;

并将线设置向量y替换为

arma::vec y = arma::conv_to<stdvec>::from(Y.row(i)); 

在我的qSelectMbyrow 函数中，虽然该函数随后运行，但与在列上应用相比，它运行缓慢，并且如果我运行 100 次，我的 R 会话也会崩溃。

【问题讨论】：

你试过把arma::vec y = Y.row(i);改成arma::rowvec y = Y.row(i);吗？

哈，谢谢-成功了！如果您愿意，可以将其发布为答案！

很高兴就这么简单——添加了答案

为了更快的速度，你可以省略arma::mat Y(M.memptr(), M.n_rows, M.n_cols);完成的矩阵复制已经有arma::rowvec y = Y.row(i);正在制作的复制

谢谢你让我知道 - 这确实也可以正常工作（也可以与 OpenMP 结合使用）并且速度稍快！

【参考方案1】：

问题在于arma::vec 实际上是arma::colvec（请参阅the docs）。所以，我们可以通过改变来解决这个问题

arma::vec y = Y.row(i);

（这是不兼容的，因为它认为你想要一个一列的矩阵，但你试图给它一个一行的矩阵）到

arma::rowvec y = Y.row(i);