r - 矩阵列表乘以标量，Rcpp 中不保留维度属性

Question

我正在使用 Rcpp 来加速一些 R 代码(实际上这是我 2014 年“待办事项”列表中的项目之一)，部分代码包括将矩阵列表乘以标量，我能够得到结果，尽管如此，矩阵不再是矩阵，而是向量，我想要一个矩阵列表作为最终输出。

这是我到目前为止的代码:

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

// I got this template from here: http://stackoverflow.com/a/18014655/1315767
template <typename WHAT>
class ListOf : public List {
public:
    template <typename T>
    ListOf( const T& x) : List(x){}

    WHAT operator[](int i){ return as<WHAT>( ( (List*)this)->operator[]( i) ) ; }

} ;

// [[Rcpp::export]]


List FooList(NumericVector fi1, ListOf<NumericMatrix> Ct){

  List TempList(Ct.size());
  NumericMatrix ct(2,2);


  for(int i=0; i<Ct.size(); i++){
    ct = Ct[i] ;
    TempList[i] = ct * fi1[i] ;  // multiply each matrix by each scalar in fi1
  }
   return TempList;
}

运行此代码时，我得到以下信息:

> sourceCpp("FooList.cpp")
> A <- replicate(3,matrix(1:4, 2), simplify=FALSE) # a list of matrices
> vec <- 0.5 * c(1:3)                              # a vector
> FooList(vec, A)  # dim are not preserved
[[1]]
[1] 0.5 1.0 1.5 2.0

[[2]]
[1] 1 2 3 4

[[3]]
[1] 1.5 3.0 4.5 6.0

FooList 给出的输出可以，但格式不是，我希望得到这样的东西:

[[1]]
     [,1] [,2]
[1,]  0.5  1.5
[2,]  1.0  2.0

[[2]]
     [,1] [,2]
[1,]    1    3
[2,]    2    4

[[3]]
     [,1] [,2]
[1,]  1.5  4.5
[2,]  3.0  6.0

我不明白为什么我会得到这个输出，因为 ct是一个矩阵，如果我去掉 fi1[i]输出确实是一个矩阵列表，我什至尝试使用 as_scalar(fi)我和以前一样。我也尝试使用 ct.attr("dim") = Dimension(2, 2);没有成功。

Answer 1

关键问题是，当您在 C++ 中将矩阵乘以标量时，您使用的是 Rcpp * 的语法糖，这是矢量化的。无论出于何种原因，它都不了解如何返回矩阵(我没有广泛查看文档)。

如果我们改为将每个矩阵的每个元素乘以标量，就会得到预期的结果:

FooList.R

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

// I got this template from here: http://stackoverflow.com/a/18014655/1315767
template <typename WHAT>
class ListOf : public List {
public:
    template <typename T>
    ListOf( const T& x) : List(x){}

    WHAT operator[](int i){ return as<WHAT>( ( (List*)this)->operator[]( i) ) ; } 

} ;

// [[Rcpp::export]]


List FooList(NumericVector fi1, ListOf<NumericMatrix> Ct){

  List TempList(Ct.size());
  NumericMatrix ct(2,2);


  for(int i=0; i<Ct.size(); i++){
    ct = Ct[i] ;
    for (int j=0; j < ct.nrow(); j++) {
      for (int k=0; k < ct.ncol(); k++) {
        ct(j, k) *= fi1[i];  // Multiply each element of the matrix by the scalar in fi1
      }   
    }   
    TempList[i] = ct;                                                                                                                                                                                                                      
  }
   return TempList;
}

互动环节:

> sourceCpp("FooList.cpp")
> A <- replicate(3,matrix(1:4, 2), simplify=FALSE) # a list of matrices
> vec <- 0.5 * c(1:3)                              # a vector
> FooList(vec, A) 
[[1]]
     [,1] [,2]
[1,]  0.5  1.5
[2,]  1.0  2.0

[[2]]
     [,1] [,2]
[1,]    1    3
[2,]    2    4

[[3]]
     [,1] [,2]
[1,]  1.5  4.5
[2,]  3.0  6.0