我stats::uniroot在data.table 中的一百萬行上運行函式。這是一個玩具示例-
library(data.table)
cumhaz <- function(t, a, b) b * (t/b)^a
froot <- function(x, u, a, b) cumhaz(x, a, b) - u
n <- 50000
u <- -log(runif(n))
a <- 1/2
b <- 1
dt = data.table(u = u, a = a, b = b)
print(system.time(
dt[, c := uniroot(froot, u=u, a=a, b=b, interval= c(0.01, 10), extendInt="yes")$root, by = u]
))
在上面的代碼中，50,000 行花費的時間接近 8 秒。
有沒有更快的替代uniroot函式可以大大減少這個時間？
uj5u.com熱心網友回復：
160 秒 (1e6/5e4 * 8) 對于一百萬行來說對我來說聽起來還不錯（盡管您的實際功能可能比froot您在這里使用的要慢得多？）。這可以簡單地并行化，在不同的內核上運行不同的塊（參見例如這個問題的答案）。
你有多需要extendInt？如果我制作一個uniroot()只有核心功能的函式的黑客版本，沒有引數測驗邏輯等，我可以將速度提高三倍。 但是，如果您的目標函式比你在這里給出的例子；如果是這種情況，您應該專注于加速您的目標函式（我嘗試froot在 C 中通過重新編碼您的代碼Rcpp，但在這種情況下它并沒有真正幫助 - 該函式非常微不足道，函式呼叫開銷大??部分時間...）
我只用 5000 行做了這個，以便于進行基準測驗：
n <- 5000
u <- -log(runif(n))
a <- 1/2
b <- 1
dt = data.table(u = u, a = a, b = b)
最小功能：
uu <- function(f, lower, upper, tol = 1e-8, maxiter =1000L, ...) {
f.lower <- f(lower, ...)
f.upper <- f(upper, ...)
val <- .External2(stats:::C_zeroin2, function(arg) f(arg, ...),
lower, upper, f.lower, f.upper, tol, as.integer(maxiter))
return(val[1])
}
檢查我們是否得到相同的結果：
identical(uniroot(froot, u = 3.242, a=0.5, b=1, interval = c(0.01,100))$root,
uu(froot, u = 3.242, a=0.5, b=1, lower = 0.01, upper = 100))
## TRUE
基準測驗包；將評估包裝在函式中以保持緊湊
library(rbenchmark)
f1 <- function() {
dt[, c := uniroot(froot_cpp, u=u, a=a, b=b, interval= c(0.01, 10), extendInt="yes")$root, by = u]
}
f2 <- function() {
dt[, c := uu(froot, u=u, a=a, b=b, lower = 0.01, upper = 100), by = u]
}
bb <- benchmark(f1(), f2(),
columns =c("test", "replications", "elapsed", "relative"))
結果：
test replications elapsed relative
1 f1() 100 34.616 3.074
2 f2() 100 11.261 1.000
uj5u.com熱心網友回復：
請注意，所示函式的逆可以顯式計算為
f2 <- function(x) (b^a * x / b)^(1/a)
a <- 1/2
b <- 1
all.equal(f(.5), f2(.5)) # f defined below using uniroot
## [1] TRUE
然而，假設實際上你有一個更復雜的函式，我們可以使用切比雪夫近似來得到它的近似值。請注意，a 和 b 是問題中的常量，因此我們也假設是下面的情況，即 f 使用在全域環境中設定的常量 a 和 b。以下代碼的運行速度比基準問題中的代碼快 100 倍，使用 9 次多項式，并且與 uniroot 給出的答案相差 1e-4。如果您需要更高的準確性，請使用更高的學位。
library(data.table)
library(pracma)
set.seed(123)
cumhaz <- function(t, a, b) b * (t/b)^a
froot <- function(x, u, a, b) cumhaz(x, a, b) - u
n <- 5000
u <- -log(runif(n))
a <- 1/2
b <- 1
dt = data.table(u = u, a = a, b = b)
dt2 <- copy(dt)
f <- function(u) {
uniroot(froot, u=u, a=a, b=b, interval= c(0.01, 10), extendInt="yes")$root
}
library(microbenchmark)
microbenchmark(times = 10,
orig = dt[, c := uniroot(froot, u=u, a=a, b=b, interval= c(0.01, 10), extendInt="yes")$root, by = u],
cheb = dt2[, c := chebApprox(u, Vectorize(f), min(u), max(u), 9)]
)
## Unit: milliseconds
## expr min lq mean median uq max neval cld
## orig 943.5323 948.9321 961.00361 958.91970 972.6308 982.0060 10 b
## cheb 9.3752 9.7513 10.67386 10.02555 10.3411 16.9475 10 a
max(abs(dt$c - dt2$c))
## [1] 8.081021e-05
uj5u.com熱心網友回復：
確切的問題有很好的答案，但有一些關于一般 R 實踐的說明。
當順序無關緊要時使用
在 OP 中，我們使用的by = u是每次運行一行。這是低效的！data.table將對您進行排序，確定分組，并且由于它們是真正的非常隨機的數字，因此最終會得到與行一樣多的分組。
相反，我們可以使用Map()或mapply()遍歷將提高性能的行。請注意，不清楚是否a以及b實際上是否因行而異 - 如果它們確實是常量，我們可能希望將它們從 data.table 中取出并將它們作為常量傳遞。
uniroot2 = function(...) uniroot(...)$root ## helper function
dt[, c2 := mapply(uniroot2, u, a,b,
MoreArgs = list (f = froot,
interval = c(0.01, 10),
extendInt = 'yes'))]
## for n = 5000
## # A tibble: 2 x 13
## expression min median `itr/sec` mem_alloc `gc/sec` n_itr n_gc total_time
## <bch:expr> <bch:t> <bch:t> <dbl> <bch:byt> <dbl> <int> <dbl> <bch:tm>
## 1 OP 1.17s 1.17s 0.851 170KB 2.55 1 3 1.17s
## 2 no_by 857.2ms 857.2ms 1.17 214KB 3.50 1 3 857.2ms
##
## Warning message:
## Some expressions had a GC in every iteration; so filtering is disabled.
請注意，一旦我們在 中設定了它mapply，就可以輕松future.apply::future_mapply()地并行化我們的呼叫。這比no_by我的筆記本電腦上的示例快 2.5 倍。
library(future.apply)
plan(multisession)
dt[, c3 := future_mapply(uniroot2, u, a,b,
MoreArgs = list (f = froot,
interval = c(0.01, 10),
extendInt = 'yes')
, future.globals = "cumhaz")] ## see next section for how we could remove this
函式呼叫需要時間
在您的示例中，您將兩個函式定義為：
cumhaz <- function(t, a, b) b * (t/b)^a
froot <- function(x, u, a, b) cumhaz(x, a, b) - u
當性能是一個問題并且簡化起來很簡單時，您可能想要簡化。
froot2 = function(x, u, a, b) b * (x / b) ^ a - u
超過一百萬個回圈，額外的呼叫cumhaz()加起來：
x = 2.5; u = 1.5; a = 0.5; b = 1
bench::mark(froot_rep = for (i in 1:1e6) {froot(x=x, u=u, a=a, b=b)},
froot2_rep = for (i in 1:1e6) {froot2(x=x, u=u, a=a, b=b)})
## # A tibble: 2 x 13
## expression min median `itr/sec` mem_alloc `gc/sec` n_itr n_gc total_time
## <bch:expr> <bch:t> <bch:t> <dbl> <bch:byt> <dbl> <int> <dbl> <bch:tm>
## 1 froot_rep 4.74s 4.74s 0.211 13.8KB 3.38 1 16 4.74s
## 2 froot2_rep 3.17s 3.17s 0.315 13.8KB 2.84 1 9 3.17s
##
## Warning message:
## Some expressions had a GC in every iteration; so filtering is disabled.
因為uniroot會進一步增加呼叫，默認最大迭代次數為 1,000！這意味著cumhaz()在優化程序中我們的成本在 1.5 到 1,500 秒之間。而作為@G。Grothendieck 指出，有時我們實際上可以直接求解和使用直接向量化的方法，而不是依賴uniroot或optimize。
轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/350432.html
上一篇：如何通過列中的近似值提取行？