RcppEigen - What's wrong with this matrix multiplication?

eigen, matrix, r, rcpp

Solution

Thanks for your edit. Multiplying a 5 x 5 matrix with 5 x 1 vector could never produce a 5 x 5. You do want the identity matrix here, which is what `diag(rep(1,5))` gives you, as does `diag(5)`.

A brief search at the Eigen docs suggests something like the following using `MatrixXd::Identity()` should do:

#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);
}

Just run `sourceCpp("nameOfTheFile.cpp")` on it, and then:

 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>

with the same `tm` as you had.

Edit 1: For the actual Hadamard product, you may just want to copy, say, the routine from Octave or find another one somewhere...

Edit 2: More concise interface placing `Map<MatrixXd>` directly in function signature.

Problem

I'm trying for a simple Hadamard product, i.e. a matrix in which the first element of vector `fv` is multiplied by all column-1 elements of matrix `tm`, the second by column-2, etc. Minimal example: ``` 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 ``` This should simply return `tm` rather than broadcasting the first column, I don't know what it actually thinks it is doing. Am I omitting something obvious? EDIT: `etest(tm,diag(fv))` gives me what I want, but this should be doable from within eigen?

Original source