logspace add RTMB

[Jan-Ole Koslik]

Quick question: Is there a way to do logspace add for more than 2 numbers with RTMB? 

I'm not sure how do have a custom implementation without using max().

Cheers

[Kasper Kristensen]

MakeTape(max, 1:10)

works for me. Maybe it's not in the CRAN version yet?

[Ben Bolker]

Stupid question, how can max(advector) possibly be differentiable ... ??

It's a little inefficient, but a loop would seem to work:

library(RTMB)
library(statnet.common)
library(microbenchmark)

set.seed(101)
xx <- rnorm(100, mean = 1, sd = 2)
r1 <- log_sum_exp(xx)

myfun <- function(x) {
res <- logspace_add(x[1], x[2])
if (length(x)>2) {
for (i in 3:length(x)) {
res <- logspace_add(res, x[i])
}
}
res
}
all.equal(myfun(xx), r1)
ff <- MakeADFun(myfun, xx)
ff$fn()

microbenchmark(ff$fn(xx), myfun(xx), log_sum_exp(xx))

Unit: microseconds
expr min lq mean median uq
ff$fn(xx) 11.521 12.1875 15.16938 15.043 16.581
myfun(xx) 2641.895 2680.3430 2759.94408 2700.575 2722.672
log_sum_exp(xx) 1.453 1.6285 2.15878 1.994 2.370

> --
> To post to this group, send email to us...@tmb-project.org. Before
> posting, please check the wiki and issuetracker at https://github.com/
> kaskr/adcomp/ <https://github.com/kaskr/adcomp/>. Please try to create a
> simple repeatable example to go with your question (e.g issues 154, 134,
> 51). Use the issuetracker to report bugs.
> ---
> You received this message because you are subscribed to the Google
> Groups "TMB Users" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to tmb-users+...@googlegroups.com <mailto:tmb-
> users+un...@googlegroups.com>.
> To view this discussion visit https://groups.google.com/d/msgid/tmb-
> users/564fbde3-e13f-48b1-887a-fa84c6e7556fn%40googlegroups.com <https://
> groups.google.com/d/msgid/tmb-users/564fbde3-e13f-48b1-887a-
> fa84c6e7556fn%40googlegroups.com?utm_medium=email&utm_source=footer>.

--
Dr. Benjamin Bolker
Professor, Mathematics & Statistics and Biology, McMaster University
Director, School of Computational Science and Engineering
> E-mail is sent at my convenience; I don't expect replies outside of
working hours.

[Kasper Kristensen]

Ben, although max isn't differentiable it's OK as long as the result is:

lse <- function(x) { m <- max(x); log(sum(exp(x-m))) + m }

This is faster to tape:

system.time(MakeTape(lse, numeric(1e6)))

system.time(MakeTape(myfun, numeric(1e6)))

[Jan-Ole Koslik]

Thanks to both of you!