Rcpparmadillo 矩阵产品性能

Question

有人可以向我解释为什么在我的代码中添加 arma::mat P(X * arma::inv(X.t() * X) * X.t()); 后计算变得如此缓慢。上次我对代码进行基准测试时，平均值增长了 164 倍。

// [[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)
  */

此致，雅各布

Answer 1

好问题。不完全确定为什么除了我所做的一些笔记之外速度会增加。所以，请注意。

考虑这里使用的 n 是 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 条目表明：

if matrix A is know to be symmetric positive definite, using inv_sympd() is faster

if matrix A is know to be diagonal, use inv( diagmat(A) )

to solve a system of linear equations, such as Z = inv(X)*Y, using solve() is faster and more accurate

在这种情况下，第三点特别有趣，因为它为 inv(X.t() * X)*X.t() => solve(X.t() * X, X.t())

提供了更优化的例程

Answer 2

你所做的非常接近 fastLm() 我多年来修改了很多次。从中我们可以得出一些结论：

不要直接X(X' X)^1 X'。使用 solve().
永远不要处理公式对象。对 X 和 y.

这里benchmark example说明了解析公式如何破坏矩阵代数的所有收益。

顺便说一句，R本身有针对秩亏矩阵的旋转操作。这有助于变形矩阵；在许多 "normal" 情况下你应该没问题。

Rcpparmadillo 矩阵产品性能

Rcpparmadillo matrixproduct performance

r

rcpp