RcppEigen - 这个矩阵乘法有什么问题？

Question

我正在尝试一个简单的 Hadamard 产品，ie 一个矩阵，其中向量 fv 的第一个元素乘以矩阵 tm 的所有第 1 列元素，第二个乘以第 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)) 给了我我想要的东西，但这应该可以在 eigen 内实现？

Answer 1

感谢您的编辑。将 5 x 5 矩阵与 5 x 1 向量相乘永远不会产生 5 x 5。您确实需要这里的单位矩阵，这就是 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>直接放在函数签名中。