2014-07-07 73 views
1

我試圖爲一個簡單的Hadamard乘積,即其中矢量fv的第一個元素是由矩陣tm的所有列-1元素相乘的矩陣,所述第二通過柱-2等RcppEigen - 這種矩陣乘法有什麼問題?

小例子:

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))給了我我想要的,但這應該是可以從特徵內實現嗎?

回答

3

感謝您的編輯。將5×5矩陣乘以5×1向量不能產生5×5。你確實需​​要這裏的單位矩陣,這就是diag(rep(1,5))給你的結果,diag(5)也是如此。

在本徵文檔的簡單搜索提出像下面使用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:對於實際阿達馬的產品,你可能只是想複製,比如,從八度的常規或另找一個地方......

編輯2:更簡潔的界面將Map<MatrixXd>直接在函數簽名中。

+0

謝謝德克。對不起,我不是更清楚,這解決了我的例子,但不是更普遍的問題。我真正想要做的就是傳入一個任意的向量作爲'fv'並且進行列式乘法運算。我承認它不是線性代數意義上的標準乘法,但是我希望看到在元素純Rcpp解決方案方面有所收穫,如果我可以從矩陣的角度來理解的話。我從這裏列出的「矩陣/矢量產品」中繪製:http://eigen.tuxfamily.org/dox/group__QuickRefPage.html –

+0

我們似乎在同一時間編輯過,因此請參閱我添加的兩行。如果您需要不同的操作,請根據需要實施。 –

+0

謝謝。事實上,這似乎沒有適當的執行,我濫用了操作員。我會嘗試玩更多的對角線。 –