There are great answers to the exact question, but a couple of notes on general R practices.
Using by when order doesn't matter
In the OP, we are using by = u
so that each row is run one at a time. This is inefficient! data.table
will order u, determine groupings, and since they are real very random numbers, end up with as many groupings as rows.
Instead, we can use Map()
or mapply()
to iterate through the rows which will improve performance. Note, it's unclear whether a
and b
actually vary by row - if they truly are constant, we would likely want to take them out of the data.table and pass them as constants.
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.
Note, once we have it set up in mapply
, it is trivial to use future.apply::future_mapply()
to parallelize our call. This is 2.5 times faster than the no_by
example above on my laptop.
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
Function calls take time
In your example, you define two functions as:
cumhaz <- function(t, a, b) b * (t/b)^a
froot <- function(x, u, a, b) cumhaz(x, a, b) - u
When performance is an issue and it is trivial to simplify, you may want to simplify.
froot2 = function(x, u, a, b) b * (x / b) ^ a - u
Over a million of loops, the additional call to cumhaz()
adds up:
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.
Since uniroot
would further increase the calls with a default max iterations of 1,000! That means cumhaz()
costs us somewhere between 1.5s and 1,500s during the optimization. And as @G. Grothendieck pointed out, sometimes we can actually directly solve and used direct vectorized methods instead of relying on uniroot
or optimize
.