Производительность матрицы Rcpparmadillo

0 Jakob Jul Elben [2017-06-22 11:58:00]

Может кто-нибудь объяснить мне, почему вычисления становятся намного медленнее, когда я добавляю arma::mat P(X * arma::inv(Xt() * X) * Xt()); к моему коду. Среднее значение выросло с коэффициентом 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)
  */

С наилучшими пожеланиями, Якоб

r rcpp

Оригинал вопроса (на английском)

1 Решение Dirk Eddelbuettel [2017-06-22 21:09:00]

То, что вы делаете, очень близко к fastLm() которое я пересматривал много раз за эти годы. Из этого можно сделать несколько выводов:

Не произвольно X (X 'X) ^ 1 X'. Используйте solve().
Никогда не работайте с объектом формулы. Используйте матрицу и вектор для X и y.

Вот примерный пример, иллюстрирующий, как разбор формулы разрушает все выгоды от матричной алгебры.

В стороне, R сама по себе поворачивает операции для ранга-дефицитной матрицы. Это помогает с деформированными матрицами; во многих "нормальных" случаях вы должны быть в порядке.

Оригинал ответа (на английском)

1 coatless [2017-06-22 18:48:00]

Отличный вопрос. Не совсем уверен, почему увеличение скорости за пределами нескольких заметок, которые я сделал. Поэтому будьте осторожны.

Рассмотрим, что 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 указано, что:

если матрица A известна как симметричная положительно определенная, с использованием inv_sympd() выполняется быстрее

если матрица А известна как диагональная, используйте inv (diagmat (A))

для решения системы линейных уравнений, таких как Z = inv (X) * Y, с использованием solve() быстрее и точнее

Третья точка особенно интересна в этом случае, поскольку она дает более оптимизированную подпрограмму для inv(Xt() * X)*Xt() => solve(Xt() * X, Xt())

Оригинал ответа (на английском)

Производительность матрицы Rcpparmadillo

2 ответа