Rcpp：* 的语法糖在处理 NumericMatrix 时会产生意想不到的结果答案

【问题标题】：Rcpp: Syntactic sugar for * produces unexpected results when dealing with NumericMatrixRcpp：* 的语法糖在处理 NumericMatrix 时会产生意想不到的结果
【发布时间】：2014-01-06 13:18:46
【问题描述】：

A recently asked question 让我相信Rcpp 的* 的语法糖不能按预期工作。在链接的问题中，用户试图将矩阵乘以标量。

R 代码

这是我们试图在 Rcpp 中实现的目标，但现在是简单的 R：

> m <- matrix(0:3, 2, 2)
> m * 3
     [,1] [,2]
[1,]    0    6
[2,]    3    9

Rcpp 代码

我创建了一些最小的示例来演示上述问题以及一些意外行为。首先请注意，我一直使用List 作为返回类型，因为它不需要我提前声明适当的类型：

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]

List FooMat() {
  // Create a fill a 2x2 matrix
  NumericMatrix tmp(2,2);
  for (int i = 0; i < 4; i++) {
    tmp[i] = i;
  }

  return List::create(tmp);
}

// [[Rcpp::export]]
List FooMat2() {
  // Create a fill a 2x2 matrix
  NumericMatrix tmp(2,2);
  for (int i = 0; i < 4; i++) {
    tmp[i] = i;
  }

  NumericVector x(1);                                                                                                                                                                                                                      
  x[1] = 3;

  return List::create(tmp * x); 
}

// [[Rcpp::export]]

List FooMat3() {
  // Create a fill a 2x2 matrix
  NumericMatrix tmp(2,2);
  for (int i = 0; i < 4; i++) {
    tmp[i] = i;
  }

  NumericVector x(1);
  x[1] = 3;

  return List::create(tmp * x[1]);
}

// [[Rcpp::export]]

List FooMat4() {
  // Create a fill a 2x2 matrix
  NumericMatrix tmp(2,2);
  for (int i = 0; i < 4; i++) {
    tmp[i] = i;
  }

  return List::create(tmp * 3); 
}

现在，如果我们获取文件，我们会得到一些奇怪的行为：

# Proof that we can return a NumericMatrix in a List:
> FooMat()
[[1]]
     [,1] [,2]
[1,]    0    2
[2,]    1    3

# Multiply the whole NumericMatrix by a whole NumericVector
# whose size is 1. Unsafe behaviour?
> FooMat2()
[[1]]
[1]  0.000000e+00  3.000000e+00 1.388988e-309 2.083483e-309

# Multiply the whole NumericMatrix by the first element of
# The NumericVector. Results are correct, but `*` converts
# the answer to a NumericVector instead of a NumericMatrix
> FooMat3()
[[1]]
[1] 0 3 6 9

# Same as FooMat3() except now we just multiply the NumericMatrix
# by an integer
> FooMat4()
[[1]]
[1] 0 3 6 9

第一，Rcpp 提供的* 的语法糖似乎不能正确处理矩阵与标量的乘法。第二，乘以整数NumericVector，如FooMat2() 会导致不安全的行为。

【问题讨论】：

标签： r rcpp

【解决方案1】：

正如我在之前的答案中所说，当我需要对矩阵进行实际数学运算时，我会使用 Armadillo 对象：

R> cppFunction('arma::mat scott(arma::mat x, double z) { 
+                 return(x*z); }', 
+              depends="RcppArmadillo")
R> scott(matrix(1:4,2), 2)
     [,1] [,2]
[1,]    2    6
[2,]    4    8
R>

糖操作很好，但不完整。不过，我们肯定会打补丁。

正如我们之前所说的：rcpp-devel 是正确的支持渠道。

编辑（2016 年 10 月或 2 1/2 年后）：搜索其他内容让我回到了这里。在 Rcpp 0.12.* 系列中，一些在矩阵和向量之间进行运算时会起作用，因此基本的“矩阵时间标量”现在可以按您的预期工作：

R> cppFunction("NumericMatrix testmat(NumericMatrix m, double multme) { 
+               NumericMatrix n = m * multme; 
+               return n; }") 
R> testmat(matrix(1:4,2), 1)
     [,1] [,2]
[1,]    1    3
[2,]    2    4
R> testmat(matrix(1:4,2), 3)
     [,1] [,2]
[1,]    3    9
[2,]    6   12
R>

不过，我可能仍会使用 RcppArmadillo 进行矩阵数学运算。

【讨论】：

这澄清了很多我之前的问题出了什么问题。谢谢 Dirk，当然是 +1！
这是一个棘手的问题，使用 NumericMatrix 很诱人，因为我们似乎提供了全功能支持。实际上，我们没有——但康拉德的犰狳坚如磐石、经过充分测试、完整且得到 R. Eigen 的大力支持。也是一个公平的选择。
谢谢德克！我浏览了 Rcpp 网站并没有看到支持链接，所以我回到了这里。

【解决方案2】：

这是一个糟糕的设计决策的不幸后果，即使 Rcpp 矩阵从 Rcpp 向量派生。

我可能会在我现在维护的 Rcpp 实现中恢复这个决定：Rcpp11 和 Rcpp98。我不再认为让 Matrix 从 Vector 派生有任何好处，它会妨碍this file 末尾使用的 CRTP。

【讨论】：