r - RcppEigen-此矩阵乘法有什么问题？

我正在尝试一个简单的Hadamard乘积，即一个矩阵，其中矢量fv的第一个元素乘以矩阵tm的所有column-1元素，第二个乘以column-2，依此类推。

最小示例:

set.seed(123)
tm <- matrix(rnorm(25,2,1),nrow=5)
fv <- rep(1,5)

tm

         [,1]      [,2]     [,3]       [,4]      [,5]
[1,] 1.439524 3.7150650 3.224082 3.78691314 0.9321763
[2,] 1.769823 2.4609162 2.359814 2.49785048 1.7820251
[3,] 3.558708 0.7349388 2.400771 0.03338284 0.9739956
[4,] 2.070508 1.3131471 2.110683 2.70135590 1.2711088
[5,] 2.129288 1.5543380 1.444159 1.52720859 1.3749607


library(inline)
etest <- cxxfunction(signature(tm="NumericMatrix",
                               fv="NumericVector"),
                     plugin="RcppEigen",
                     body="
NumericVector fvv(fv);
NumericMatrix tmm(tm);

const Eigen::Map<Eigen::MatrixXd> ttm(as<Eigen::Map<Eigen::MatrixXd> >(tmm));
const Eigen::Map<Eigen::VectorXd> ffv(as<Eigen::Map<Eigen::VectorXd> >(fvv));

Eigen::MatrixXd prod = ttm*ffv.transpose();
return(wrap(prod));
                     ")

etest(tm,fv)

         [,1]     [,2]     [,3]     [,4]     [,5]
[1,] 1.439524 1.439524 1.439524 1.439524 1.439524
[2,] 1.769823 1.769823 1.769823 1.769823 1.769823
[3,] 3.558708 3.558708 3.558708 3.558708 3.558708
[4,] 2.070508 2.070508 2.070508 2.070508 2.070508
[5,] 2.129288 2.129288 2.129288 2.129288 2.129288

这应该只是返回tm而不是广播第一列，我不知道它实际上在做什么。我是否遗漏了明显的内容？

编辑:etest(tm,diag(fv))给了我我想要的东西，但这应该可以在本征内实现？

最佳答案

感谢您的修改。将5 x 5矩阵与5 x 1向量相乘永远不会产生5 x5。您确实想要此处的单位矩阵，即diag(rep(1,5))和diag(5)所提供的。

在Eigen文档上进行的简短搜索建议使用MatrixXd::Identity()进行以下操作:

#include <RcppEigen.h>

using namespace Eigen;
using namespace Rcpp;

// [[Rcpp::depends(RcppEigen)]]

// [[Rcpp::export]]
MatrixXd etest(const Map<MatrixXd> ttm));
  const MatrixXd id = MatrixXd::Identity(ttm.rows(), ttm.cols());

  Rcout << "ttm\n " << ttm << std::endl;
  Rcout << "id\n " << id << std::endl;
  MatrixXd res = ttm*id;
  Rcout << "res\n " << res << std::endl;

  return(res);
}

只需对其运行sourceCpp("nameOfTheFile.cpp")，然后:

 R> sourceCpp("/tmp/hada.cpp")
 R> etest(tm)
 ttm
    1.43952   3.71506   3.22408   3.78691  0.932176
   1.76982   2.46092   2.35981   2.49785   1.78203
   3.55871  0.734939   2.40077 0.0333828  0.973996
   2.07051   1.31315   2.11068   2.70136   1.27111
   2.12929   1.55434   1.44416   1.52721   1.37496
 id
  1 0 0 0 0
 0 1 0 0 0
 0 0 1 0 0
 0 0 0 1 0
 0 0 0 0 1
 res
    1.43952   3.71506   3.22408   3.78691  0.932176
   1.76982   2.46092   2.35981   2.49785   1.78203
   3.55871  0.734939   2.40077 0.0333828  0.973996
   2.07051   1.31315   2.11068   2.70136   1.27111
   2.12929   1.55434   1.44416   1.52721   1.37496
         [,1]     [,2]    [,3]      [,4]     [,5]
 [1,] 1.43952 3.715065 3.22408 3.7869131 0.932176
 [2,] 1.76982 2.460916 2.35981 2.4978505 1.782025
 [3,] 3.55871 0.734939 2.40077 0.0333828 0.973996
 [4,] 2.07051 1.313147 2.11068 2.7013559 1.271109
 [5,] 2.12929 1.554338 1.44416 1.5272086 1.374961
 R>

使用与您相同的tm。

编辑1:对于实际的Hadamard产品，您可能只想从Octave复制例程，或者在某个地方找到另一个例程。

编辑2:更简洁的界面，直接将Map<MatrixXd>放置在函数签名中。