Rcpparmadillo矩陣產品性能

有人可以向我解釋爲什麼當我將arma::mat P(X * arma::inv(X.t() * X) * X.t());添加到我的代碼時，計算變得如此之慢。上次我對代碼進行基準測試時，平均數增長了164倍。Rcpparmadillo矩陣產品性能

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

#include <RcppArmadillo.h> 
using namespace Rcpp; 

//[[Rcpp::export]] 
List test1(DataFrame data, Language formula, String y_name) { 
    Function model_matrix("model.matrix"); 
    NumericMatrix x_rcpp = model_matrix(formula, data); 
    NumericVector y_rcpp = data[y_name]; 
    arma::mat X(x_rcpp.begin(), x_rcpp.nrow(), x_rcpp.ncol()); 
    arma::colvec Y(y_rcpp.begin(), y_rcpp.size()); 

    arma::colvec coef = inv(X.t() * X) * X.t() * Y; 
    arma::colvec resid = Y - X * coef; 
    arma::colvec fitted = X * coef; 

    DataFrame data_res = DataFrame::create(_["Resid"] = resid, 
        _["Fitted"] = fitted); 

    return List::create(_["Results"] = coef, 
         _["Data"] = data_res); 
} 

//[[Rcpp::export]] 
List test2(DataFrame data, Language formula, String y_name) { 
    Function model_matrix("model.matrix"); 
    NumericMatrix x_rcpp = model_matrix(formula, data); 
    NumericVector y_rcpp = data[y_name]; 
    arma::mat X(x_rcpp.begin(), x_rcpp.nrow(), x_rcpp.ncol()); 
    arma::colvec Y(y_rcpp.begin(), y_rcpp.size()); 

    arma::colvec coef = inv(X.t() * X) * X.t() * Y; 
    arma::colvec resid = Y - X * coef; 
    arma::colvec fitted = X * coef; 

    arma::mat P(X * arma::inv(X.t() * X) * X.t()); 

    DataFrame data_res = DataFrame::create(_["Resid"] = resid, 
             _["Fitted"] = fitted); 

    return List::create(_["Results"] = coef, 
         _["Data"] = data_res); 
} 

/*** R 
data <- data.frame(Y = rnorm(10000), X1 = rnorm(10000), X2 = rnorm(10000), X3 = rnorm(10000)) 
microbenchmark::microbenchmark(test1(data, Y~X1+X2+X3, "Y"), 
           test2(data, Y~X1+X2+X3, "Y"), times = 10) 
    */

最好的問候，雅各布

來源

2017-06-22 Jakob Jul Elben

你在做什麼是非常接近fastLm()這是我多年來多次修改。由此我們可以得出幾點結論：

直接不要X（X'X）^ 1 X'。使用solve()。
不要有史以來工作了一個公式對象。使用矩陣和向量爲X和y。

這裏是benchmark example這裏是benchmark example說明如何解析公式破壞矩陣代數的所有收益。另一方面，R本身對等級不足矩陣進行了操作。這有助於變形矩陣;在許多「正常」情況下，你應該沒問題。

來源

2017-06-22 18:09:39

我之所以問，是因爲我目前正在研究一個新的實現，它允許管道和靈感來自於tidyverse包。我是fastlm（）函數的忠實粉絲，但這個實現是針對那些沒有任何編程經驗的經濟學生。因此，即使由於公式效率低下而導致執行速度較慢，我想使用公式。所以你會推薦使用'''model.matrix（）'''然後使用'''fastlm（）'''？ –

我從基準測試（鏈接二）得出的結論是，如果你堅持使用一個公式，那麼你不需要爲C++方面而煩惱。 –

謝謝您的全面解答。 –

大問題。不完全確定爲什麼速度會在我做出的幾個音符之外增加。所以，被警告。

考慮在這裏所使用的ñ爲10000與p是3

讓我們來看看請求的操作。我們將與coef或beta_hat的工作開始：

Beta_[p x 1] = (X^T_[p x n] * X_[n x p])^(-1) * X^T_[p x n] * Y_[n x 1]

望着P或投影/帽子矩陣：

P_[n x n] = X_[n x p] * (X^T_[p x n] * X_[n x p])^(-1) * X^T_[p x n]

所以，這裏的N矩陣是比現有矩陣充分大。矩陣乘法通常由O（n^3）（天真的教科書乘法）控制。所以，這可能會解釋大量的時間增量。

以外的是，有涉及

(X^T_[p x n] * X_[n x p])^(-1) * X^T_[p x n]

內test2導致它被重新計算的重複計算。這裏的主要問題是反向成本最高的操作。

另外，關於使用的inv API入口表明：

如果矩陣A知道是對稱正定的，使用inv_sympd（）是更快

如果矩陣A是知道爲對角的，使用INV（diagmat（A））

求解線性方程組，如Z = INV（X）的一個系統* Y，使用解決（）是更快和更準確

第三點是特別的在這種情況下利率，因爲它給出了更優化的例程inv(X.t() * X)*X.t() =>solve(X.t() * X, X.t())

來源

2017-06-22 15:48:10 coatless

Rcpparmadillo矩陣產品性能

回答

相關問題