测量 R 代码运行时间的 5 种方法

简介

平常在撰写论文时，会需要比较算法之间的计算时间。本篇文章给出几种测量 R 代码运行时间的方法。本文是小编学习过程中的笔记，主要参考博客1，2。

1. 使用 Sys.time()

小编通常使用 Sys.time() 函数来计算时间。首先记录当前运行时刻，之后运行代码，并记录运行完成时刻的时间。最后将两者时间相减即可。下面给出一个例子：

myfunction <- function() { 
Sys.sleep(10) #系统休息 
}

start_time <- Sys.time() # 记录初始时间
myfunction() # 运行你的代码
end_time <- Sys.time() # 记录终止时间

end_time - start_time # 计算时间差
# Time difference of 10.00363 secs

2. 使用 toctoc 包

函数 tic() 和 toc() 的使用方式与前面演示的 Sys.time() 类似。但是 toctoc 包更加方便。下面给出例子：

读者可以从 GitHub 或者 CRAN 上安装：

# CRAN
install.packages("tictoc")
# GitHub
library(devtools); devtools::install_github("collectivemedia/tictoc")

计算单个代码块

library(tictoc)
tic("sleeping")
print("falling asleep...")
sleep_for_a_minute()
print("...waking up")
toc()
# [1] "falling asleep..."
# [1] "...waking up"
# sleeping: 11.076 sec elapsed

嵌套多个计时器，注意这里的第一个 tic() 对应最后一个 toc()。

tic("total")  #对应最后一个toc()
tic("data generation")
X <- matrix(rnorm(50000*1000), 50000, 1000)
b <- sample(1:1000, 1000)
y <- runif(1) + X %*% b + rnorm(50000)
toc()
tic("model fitting")
model <- lm(y ~ X)
toc()
toc()
# data generation: 3.383 sec elapsed
# model fitting: 42.452 sec elapsed
# total: 46.428 sec elapsed

3. 使用 system.time()

这个使用方式类似 Sys.time()，我们继续以 myfunction() 作为例子：

system.time({ myfunction() })
#  user    system  elapsed 
#  0.039   0.044   10.001

4. 使用 rbenchmark 包

该包中的 microbenchmark() 函数可以进行精确的测量和比较 R 表达式的执行时间。它提供了更准确的结果，可以替代 system.time() 使用。

可以通过 CRAN 或者 GitHub 安装：

# CRAN
install.packages("rbenchmark")
# GitHub
devtools::install_github("eddelbuettel/rbenchmark")

这里展示R-bloggers中的例子：使用三个计算方法比较计算线性回归系数所需的时间：

library(rbenchmark)

benchmark("lm" = {
            X <- matrix(rnorm(1000), 100, 10)
            y <- X %*% sample(1:10, 10) + rnorm(100)
            b <- lm(y ~ X + 0)$coef
          },
          "pseudoinverse" = {
            X <- matrix(rnorm(1000), 100, 10)
            y <- X %*% sample(1:10, 10) + rnorm(100)
            b <- solve(t(X) %*% X) %*% t(X) %*% y
          },
          "linear system" = {
            X <- matrix(rnorm(1000), 100, 10)
            y <- X %*% sample(1:10, 10) + rnorm(100)
            b <- solve(t(X) %*% X, t(X) %*% y)
          },
          replications = 1000,
          columns = c("test", "replications", "elapsed",
                      "relative", "user.self", "sys.self"))

#           test   replications elapsed relative user.self sys.self
# 3 linear system         1000   0.064    1.000     0.062    0.001
# 1            lm         1000   0.281    4.391     0.267    0.011
# 2 pseudoinverse         1000   0.082    1.281     0.079    0.002

观察 elapsed 的结果，可以看到 lm 方法计算时间是最长的。

注意：读者想使用该包进行计算时间，需要按照文章的形式，将自己的代码填入即可。

5. microbenchmark 包

该包是 rbenchmark 包的升级版本，不仅可以给出上面的结果，还能通过可视化的方式展示结果。

可以通过 CRAN 或者 GitHub 安装：

# CRAN
install.packages("microbenchmark")
# GitHub
remotes::install_github("joshuaulrich/microbenchmark")

该函数能够使用用户指定函数自动检查基准表达式的结果。我们再次比较三种计算线性模型系数向量的方法（重复运行了100次）。

library(microbenchmark)

set.seed(2017)
n <- 10000
p <- 100
X <- matrix(rnorm(n*p), n, p)
y <- X %*% rnorm(p) + rnorm(100)

check_for_equal_coefs <- function(values) {
  tol <- 1e-12  # 添加基准
  max_error <- max(c(abs(values[[1]] - values[[2]]),
                     abs(values[[2]] - values[[3]]),
                     abs(values[[1]] - values[[3]])))
  max_error < tol
}

mbm <- microbenchmark("lm" = { b <- lm(y ~ X + 0)$coef },
               "pseudoinverse" = {
                 b <- solve(t(X) %*% X) %*% t(X) %*% y
               },
               "linear system" = {
                 b <- solve(t(X) %*% X, t(X) %*% y)
               },
               times = 100, #重复运行了100次
               check = check_for_equal_coefs) # 添加自通检查的标准

mbm
# Unit: milliseconds
#           expr      min        lq      mean    median        uq      max neval cld
#             lm 96.12717 124.43298 150.72674 135.12729 188.32154 236.4910   100   c
#  pseudoinverse 26.61816  28.81151  53.32246  30.69587  80.61303 145.0489   100  b
#  linear system 16.70331  18.58778  35.14599  19.48467  22.69537 138.6660   100 a

注意：我们使用函数参数检查来检查三种方法返回的结果的相等性（最大误差为1e-12）。如果结果不相等，则 microbenchmark 将返回错误消息。

值得一提的是，该包还整合了 ggplot2 来绘制以上结果。你只需要使用 autoplot() 即可实现：

library(ggplot2)
autoplot(mbm)

小编有话说

本推文给出了 5 种计算代码时间的方法。如果只是简单的使用，通常小编采用第一种方法。如果需要比较不同结果，并绘图。最后一种方法或许是更好的选择。

科研相关问题可见：