Printing of "Retaping: Nold=... Nnew=..." when using expAv()

[Jan-Ole Koslik]

Hi, 

I wrote a forward algorithm for a continuous-time HMM with a potentially large number of states but a sparse generator matrix Q. 

It looks like

forward_cont <- function(pi, # initial distribution (vector)
                         Q, # sparse generator matrix
                         allprobs, # matrix of state-dependent probabilities/ densities
                         timediff # time differences between observations (vector)
                         ){
  if(is.null(dim(pi))){ # if vector, reshape to matrix -> elementwise multiplication below
    pi <- matrix(pi, nrow = 1)
  }
 
  n <- nrow(allprobs) # number of observations
 
  # initialising forward algorithm
  foo <- pi * allprobs[1, , drop = FALSE]
  sumfoo <- sum(foo)
  l <- log(sumfoo)
  foo <- foo / sumfoo

  # loop over remaining time points
  for(t in 2:n){

    # foo <- foo %*% expm(Q * timediff[t-1]) # slow naive version
    foo <- t(expAv(Q * timediff[t-1], t(foo), transpose = TRUE)) # fast version
    foo <- foo * allprobs[t, , drop = FALSE]
    sumfoo <- sum(foo)
    l <- l + log(sumfoo)
    foo <- foo / sumfoo
  }
  l
}

What I commented out is the naive version, doing v %*% expm(A), which is not efficient due to expm() not exploiting sparsity. Hence, I replaced it with t(expAv(A, t(v), transpose = TRUE)), which yields the same result but is much faster.

When including this in a likelihood function, optimisation seems to work, but in between prints of "outer mgc: ..." I get loads of prints like "Retaping: Nold=53 Nnew=52" which does not happen when using expm(). 

So my question is, what does that mean, and is it something to worry about?

Cheers

[James Thorson]

I understand that it is using "uniformization," i.e., using the current parameter values defining `A` to adaptively identify the maximum number of steps (and resulting sparsity) needed to approximate v %*% expm(A) to a desired accuracy specified as argument `tol` (based on the CDF of a Poisson distribution).  So as it updates fixed effects it changes the number of steps from `Nold` to `Nnew`, up to a maximum of `Nmax`.

--
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/. 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.
To view this discussion visit https://groups.google.com/d/msgid/tmb-users/68172d77-2b77-407d-b5d3-7918052b5210n%40googlegroups.com.

[Kasper Kristensen]

Also note that this output can be disabled using the `trace` argument - see ?expAv.