-1

设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;
}

为什么这么慢？我怎样才能加快速度？

4

1 回答 1

3

欢迎来到 Stack Overflow manju。对于未来的问题，请注意，提供一个最小的可重现示例，实际上提供该示例最符合您的利益；它帮助别人帮助你。下面是一个示例，说明如何为其他人提供示例数据以供使用：

## 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)"))

在此处输入图像描述

为什么？

正如Dirk Eddelbuettel在评论中有用地指出的那样，最终 R 和犰狳都会调用 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)

在此处输入图像描述

于 2019-11-04T14:43:19.520 回答