在RcppEigen中使用cross()方法时出错

Question

我正在尝试实现一个简单的函数，我可以从 R 调用该函数来计算一对 3D 向量的叉积。明确地说，我的意思是这个叉积。

我正在尝试通过 RcppEigen 使用 Eigen 库中的

cross()

方法。我失败的尝试目前看起来像这样：

#include <Rcpp.h>
#include <RcppParallel.h>
#include <RcppEigen.h>
#include <Eigen/Eigen>
#include <Eigen/Geometry>
using namespace Rcpp;
using namespace RcppParallel;
using Eigen::Map;
using Eigen::VectorXd;
using Eigen::Vector3d;
using Eigen::MatrixXd;
// [[Rcpp::plugins("cpp11")]]
// [[Rcpp::export]]
NumericVector vectorCrossProduct3D(NumericVector x, NumericVector y) {
    const Map<VectorXd> xEigen(as<Map<VectorXd> >(x));
    const Map<VectorXd> yEigen(as<Map<VectorXd> >(y));
    VectorXd crossProduct = xEigen.cross(yEigen);
    NumericVector crossProductR = wrap(xEigen);
    return(crossProductR);
}

当前尝试编译时会导致以下类型的错误：

error: static_assert failed due to requirement 'Matrix<double, -1, 1, 0, -1, 1>::IsVectorAtCompileTime && Matrix<double, -1, 1, 0, -1, 1>::SizeAtCompileTime == 4' "THIS_METHOD_IS_ONLY_FOR_VECTORS_OF_A_SPECIFIC_SIZE"
      EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Matrix, 4)

我还尝试将

cross()

方法的输出声明为

MatrixXd

，但这也失败了。我也尝试拨打

cross3()

，但也失败了。

我觉得我在这里错过了一些非常基本的东西，但我是 RcppEigen 的新手，无法弄清楚......我将非常感谢任何对正在发生的事情的见解

更新

因此，它似乎可以通过声明

Vector3d

类型的新变量并将

cross

应用于这些变量来工作，如下所示：

// [[Rcpp::plugins("cpp11")]]
// [[Rcpp::export]]
NumericVector vectorCrossProduct3D(NumericVector x, NumericVector y) {
    Map<VectorXd> xEigen(as<Map<VectorXd> >(x));
    Map<VectorXd> yEigen(as<Map<VectorXd> >(y));
    Vector3d xEigen3 = xEigen;
    Vector3d yEigen3 = yEigen;
    Vector3d crossProduct = xEigen3.cross(yEigen3);
    NumericVector crossProductR = wrap(crossProduct);
    return(crossProductR);
}

我猜编译器会阻止将

cross

应用于未知大小的向量。然而，当我将它与普通 R 中的叉积实现进行基准测试时，它大约慢了一倍，这似乎不太正确：

vectorCrossProduct3D_R <- function(u, v) {
    w <- c(u[2]*v[3] - u[3]*v[2],
           u[3]*v[1] - u[1]*v[3],
           u[1]*v[2] - u[2]*v[1])
    return(w)
}

library(microbenchmark)

microbenchmark(vectorCrossProduct3D(1:3, 4:6), vectorCrossProduct3D_R(1:3, 4:6))

Answer 1

我不清楚你为什么使用 Eigen。您的 R 代码表明这是一个简单的算术计算，您可以简单地将其转换为 C++。您的 R 版本本质上与

pracma::cross

执行相同的操作，但没有生产代码中所需的类型和长度检查。

vectorCrossProduct_R <- function(u, v) {
  w <- c(u[2]*v[3] - u[3]*v[2],
         u[3]*v[1] - u[1]*v[3],
         u[1]*v[2] - u[2]*v[1])
  return(w)
}

做同样的事情：

Rcpp::cppFunction('
NumericVector vectorCrossProduct_cpp(NumericVector u, NumericVector v) {
  if(v.length() != 3 || u.length() != 3) {
    stop("Input vectors must be length 3");
  }
  NumericVector out = {u[1]*v[2] - u[2]*v[1],
                       u[2]*v[0] - u[0]*v[2],
                       u[0]*v[1] - u[1]*v[0]};
  return out;
}
')

测试，我们有：

vectorCrossProduct_cpp(1:3, 4:6)
#> [1] -3  6 -3
vectorCrossProduct_R(1:3, 4:6)
#> [1] -3  6 -3

但是，基准测试将显示这两种方法之间几乎没有有意义的差异。当迭代大型向量或列表时，R 相对较慢，除非迭代在底层 C 代码中进行了向量化。通过 3 个简单计算创建单个长度为三的向量将会很快。

如果现有数据结构中存储有大量 3 向量，那么您可能会通过在 C++ 中迭代它而获得显着的好处，但对于偶尔的调用，R 版本就可以了。