"Initial value in 'vmmin' is not finite" in register_fpca

[bauer-alex]

When applying the joint registration and Gaussian FPCA scheme (register_fpca) to the derivation curves of the Berkeley child growth data, the constrOptim() optimization sometimes runs into the error "Initial value in 'vmmin' is not finite". Here a code snippet to reproduce the error, and below some more insight into the problem:

library(registr)
library(dplyr)
library(ggplot2)
library(mgcv)    # for smoothing the raw curves

# Berkeley child growth data
dat <- fda::growth

# calculate the first derivative of the (slightly smoothed) curves --------
# transform to long dataset
dat_long <- data.frame(id    = rep(c(colnames(dat$hgtm), colnames(dat$hgtf)),
                                   each = length(dat$age)),
                       index = rep(dat$age, times = ncol(dat$hgtm) + ncol(dat$hgtf)),
                       value = c(as.vector(dat$hgtm), as.vector(dat$hgtf)))

# slightly smooth the raw curves
smooth_list <- lapply(unique(dat_long$id), function(curve_id) {
  d <- dat_long %>% filter(id == curve_id)
  m <- gam(value ~ s(index, bs = "cr", k = 15), data = d)
  d$value <- unname(predict.gam(m))
  return(d)
})
dat_smooth <- dplyr::bind_rows(smooth_list)

# take the first derivative
deriv_list <- lapply(unique(dat_long$id), function(curve_id) {
  d <- dat_smooth %>% filter(id == curve_id)
  data.frame(id               = curve_id,
             index            = d$index[2:nrow(d)],
             value            = diff(d$value) / diff(d$index),
             stringsAsFactors = FALSE)
})
dat_deriv <- dplyr::bind_rows(deriv_list) %>%
  mutate(id = factor(id))

# plot the first derivative
ggplot(dat_deriv, aes(x = index, y = value, col = id)) +
  geom_line() +
  ggtitle("First derivative (based on smoothed raw curves)") +
  theme(legend.position = "none")

# register the first derivative curves ------------------------------------
res <- register_fpca(Y      = dat_deriv,
                     family = "gaussian",
                     npc    = 4,
                     Kt     = 8,
                     Kh     = 4,
                     max_iterations = 100)
# -> error in iteration 13: 'Initial value in "vmmin" is not finite'

The error occurred in the joint iteration 13 for curve 12, i.e. id boy12. The curve of this id looks like this:

dat_deriv %>%
  filter(id == "boy12") %>%
  ggplot(aes(x = index, y = value)) +
  geom_line() +
  ylim(range(dat_deriv$value))

What I found out so far:

• the problem only occurrs in instances where beta_i has two nearly identical neighboring values. For curve 12 in iteration 13 beta_i equals c(7.752823, 7.752823), with beta_i[2] - beta_i[1] = 8.881784e-16.
• for curve 12, iteration 13, the corresponding gradient for the next step seems quite reasonable with values c(1.224627, 1.507873)
• the problem might only occur for really nearly identical beta_i values. E.g., in iteration 13, curve 9 has the starting values c(6.089687, 6.089687) with beta_i[2] - beta_i[1] = 4.179768e-11, causing no error in the registration step.
• the vmmin error is specific to method = "BFGS" in constrOptim(). In case it helps, here are the source codes for and (even if these source code versions already date some years back).
• as a more random side note: For iteration 13, curve 12, the problem does not occur when using constrOptim(method = "CG"). However, a similar not finite error then occurs for other curves when optimizing with CG.

[bauer-alex]

Additional notes:

• The above example only leads to an error with the current GitHub version of registr. The CRAN version 1.0.0 converges after the 12th iteration and does not throw an error.
• Also, the error only occurs with gradient = TRUE.

I will take a deeper look into this (and the penalize_beta branch, as recommended by @julia-wrobel) after the summer, probably in October.

[bauer-alex]

I could resolve all issues regarding optimization problems like the above one. These problems arose since in the joint iteration the parameter vector passed to constrOptim was sometimes slightly outside the allowed range, or slightly improper / nonmonotone as of numerical problems of constrOptim.

In branch issue_11 I now included a new function ensure_proper_beta that simply slightly alters the parameter vector to resolve these issues before the optimization call. Since the alterations are really minimal this doesn't change the results.

The branch is based on the gfpca_twoStep branch. So we should first integrate the PR for the latter branch, and then create a PR for this one.

[julia-wrobel]

Thanks Alex- gfpca_twStep has been merged with the master branch. Is the issue_11 branch ready to be merged? If so I'll run some local checks and pull it into master.

[bauer-alex]

Yes, this branch is ready to be merged.

[julia-wrobel]

Thanks Alex- it's been merged with the master branch.