RcppArmadillo: умножение диагональной матрицы очень медленно - PullRequest
Новая
       96

RcppArmadillo: умножение диагональной матрицы очень медленно

0 голосов
/ 03 ноября 2019

Пусть x - вектор, а M - матрица.

В R я могу сделать 

D <- diag(exp(x))
crossprod(M, D%M)

, а в RcppArmadillo у меня есть следующее, которое намного медленнее. 

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
arma::mat multiple_mnv(const arma::vec& x, const arma::mat& M) {
  arma::colvec diagonal(x.size())
  for (int i = 0; i < x.size(); i++)
  {
    diagonal(i) = exp(x[i]);
  }
  arma::mat D = diagmat(diagonal);
  return     M.t()*D*M;
}

Почему это так медленно? Как я могу ускорить это?

1 Ответ

3 голосов
/ 04 ноября 2019

Добро пожаловать в стек переполнения Манджу . При возникновении вопросов, пожалуйста, имейте в виду, что ожидается минимальный воспроизводимый пример , который действительно в ваших интересах;это помогает другим помочь вам. Вот пример того, как вы могли бы предоставить пример данных для работы с другими: 

## Set seed for reproducibility
set.seed(123)
## Generate data
x <- rnorm(10)
M <- matrix(rnorm(100), nrow = 10, ncol = 10)
## Output code for others to copy your objects
dput(x)
dput(M)

Это данные, с которыми я буду работать, чтобы показать, что ваш код C ++ на самом деле не медленнее чем R. Я использовал ваш код C ++ (добавив пропущенную точку с запятой): 

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
arma::mat foo(const arma::vec& x, const arma::mat& M) {
    arma::colvec diagonal(x.size());
    for ( int i = 0; i < x.size(); i++ )
    {
        diagonal(i) = exp(x[i]);
    }
    arma::mat D = diagmat(diagonal);
    return M.t() * D * M;
}

Обратите внимание, что мне пришлось сделать некоторые из моих собственных выборов относительно типа возвращаемого объекта и типов аргументов функции(это одно из мест, где вам может помочь минимальный воспроизводимый пример: что, если эти варианты повлияют на мои результаты?) Затем я создаю функцию R, которая выполняет действия foo(): 

bar <- function(v, M) {
    D <- diag(exp(v))
    return(crossprod(M, D %*% M))
}

Обратите внимание такжечто я должен был исправить твою опечатку, изменив D%M на D %*% M. Давайте еще раз проверим, дают ли они одинаковые результаты: 

all.equal(foo(x, M), bar(x, M))
# [1] TRUE

Теперь давайте посмотрим, насколько они быстры: 

library(microbenchmark)
bench <- microbenchmark(cpp = foo(x, M), R = foo(x, M), times = 1e5)
bench
# Unit: microseconds
#  expr    min     lq     mean median     uq      max
#   cpp 22.185 23.015 27.00436 23.204 23.461 31143.30
#     R 22.126 23.028 25.48256 23.216 23.475 29628.86

Они выглядят примерно одинаково для меня! Мы также можем взглянуть на график плотности времени (отбрасывая экстремальные значения, чтобы сделать вещи немного яснее): 

cpp_times <- with(bench, time[expr == "cpp"])
R_times   <- with(bench, time[expr == "R"])
cpp_time_dens <- density(cpp_times[cpp_times < quantile(cpp_times, 0.95)])
R_time_dens   <- density(R_times[R_times < quantile(R_times, 0.95)])
plot(cpp_time_dens, col = "blue", xlab = "Time (in nanoseconds)", ylab = "",
     main = "Comparing C++ and R execution time")
lines(R_time_dens, col = "red")
legend("topright", col = c("blue", "red"), bty = "n", lty = 1,
       legend = c("C++ function (foo)", "R function (bar)"))

enter image description here

Почему?

Как подсказывает Dirk Eddelbuettel в комментариях, в конце концов и R, и Armadillo все равно будут вызывать подпрограмму LAPACK или BLAS - вы не должны ожидатьбольшая разница, если вы не сможете дать Армадилло подсказку о том, как быть более эффективным.

Можем ли мы сделать код Армадилло быстрее?

Да! Как указано в комментариях mtall , мы можем дать Армадилло намек на то, что мы имеем дело с диагональной матрицей. Давай попробуем;мы будем использовать следующий код: 

// [[Rcpp::export]]
arma::mat baz(const arma::vec& x, const arma::mat& M) {
    return M.t() * diagmat(arma::exp(x)) * M;
}

и отметим его: 

all.equal(foo(x, M), baz(x, M))
# [1] TRUE
library(microbenchmark)
bench <- microbenchmark(cpp = foo(x, M), R = foo(x, M),
                        cpp2 = baz(x, M), times = 1e5)
bench
# Unit: microseconds
#  expr    min     lq     mean median     uq      max
#   cpp 22.822 23.757 27.57015 24.118 24.632 26600.48
#     R 22.855 23.771 26.44725 24.124 24.638 30619.09
#  cpp2 20.035 21.218 25.49863 21.587 22.123 36745.72

Мы видим небольшое, но верное улучшение;давайте посмотрим графически, как мы делали раньше: 

cpp_times  <- with(bench, time[expr == "cpp"])
cpp2_times <- with(bench, time[expr == "cpp2"])
R_times    <- with(bench, time[expr == "R"])
cpp_time_dens  <- density(cpp_times[cpp_times < quantile(cpp_times, 0.95)])
cpp2_time_dens <- density(cpp2_times[cpp2_times < quantile(cpp2_times, 0.95)])
R_time_dens    <- density(R_times[R_times < quantile(R_times, 0.95)])
xlims <- range(c(cpp_time_dens$x, cpp2_time_dens$x, R_time_dens$x))
ylims <- range(c(cpp_time_dens$y, cpp2_time_dens$y, R_time_dens$y))
ylims <- ylims * c(1, 1.15)
cols  <- c("#0072b2", "#f0e442", "#d55e00")
cols  <- c("#e69f00", "#56b4e9", "#009e73")
labs  <- c("C++ original", "C++ improved", "R")
plot(cpp_time_dens, col = cols[1], xlim = xlims, ylim = ylims,
     xlab = "Time (in nanoseconds)", ylab = "",
     main = "Comparing C++ and R execution time")
lines(cpp2_time_dens, col = cols[2])
lines(R_time_dens, col = cols[3])
legend("topleft", col = cols, bty = "n", lty = 1, legend = labs, horiz = TRUE)

enter image description here

...