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矩阵产品性能

回答

相关问题