Non-linear fitting with nls() is giving me singular gradient matrix at initial parameter estimates. Why?

This is my first attempt at fitting a non-linear model in R, so please bear with me.

Problem

I am trying to understand why nls() is giving me this error:

Error in nlsModel(formula, mf, start, wts): singular gradient matrix at initial parameter estimates

Hypotheses

From what I've read from other questions here at SO it could either be because:

my model is discontinuous, or
my model is over-determined, or
bad choice of starting parameter values

So I am calling for help on how to overcome this error. Can I change the model and still use nls(), or do I need to use nls.lm from the minpack.lm package, as I have read elsewhere?

My approach

Here are some details about the model:

the model is a discontinuous function, a kind of staircase type of function (see plot below)
in general, the number of steps in the model can be variable yet they are fixed for a specific fitting event

MWE that shows the problem

Brief explanation of the MWE code

step_fn(x, min = 0, max = 1): function that returns 1 within the interval (min, max] and 0 otherwise; sorry about the name, I realize now it is not really a step function... interval_fn() would be more appropriate I guess.
staircase(x, dx, dy): a summation of step_fn() functions. dx is a vector of widths for the steps, i.e. max - min, and dy is the increment in y for each step.
staircase_formula(n = 1L): generates a formula object that represents the model modeled by the function staircase() (to be used with the nls() function).
please do note that I use the purrr and glue packages in the example below.

Code

step_fn <- function(x, min = 0, max = 1) {

  y <- x
  y[x > min & x <= max] <- 1
  y[x <= min] <- 0
  y[x > max] <- 0

  return(y)
}

staircase <- function(x, dx, dy) {

  max <- cumsum(dx)
  min <- c(0, max[1:(length(dx)-1)])
  step <- cumsum(dy)

  purrr::reduce(purrr::pmap(list(min, max, step), ~ ..3 * step_fn(x, min = ..1, max = ..2)), `+`)
}


staircase_formula <- function(n = 1L) {

  i <- seq_len(n)
  dx <- sprintf("dx%d", i)

  min <-
    c('0', purrr::accumulate(dx[-n], .f = ~ paste(.x, .y, sep = " + ")))
  max <- purrr::accumulate(dx, .f = ~ paste(.x, .y, sep = " + "))

  lhs <- "y"
  rhs <-
    paste(glue::glue('dy{i} * step_fn(x, min = {min}, max = {max})'),
          collapse  = " + ")

  sc_form <- as.formula(glue::glue("{lhs} ~ {rhs}")) 

  return(sc_form)
}


x <- seq(0, 10, by = 0.01)
y <- staircase(x, c(1,2,2,5), c(2,5,2,1)) + rnorm(length(x), mean = 0, sd = 0.2)

plot(x = x, y = y)
lines(x = x, y = staircase(x, dx = c(1,2,2,5), dy = c(2,5,2,1)), col="red")


my_data <- data.frame(x = x, y = y)
my_model <- staircase_formula(4)
params <- list(dx1 = 1, dx2 = 2, dx3 = 2, dx4 = 5,
               dy1 = 2, dy2 = 5, dy3 = 2, dy4 = 1)

m <- nls(formula = my_model, start = params, data = my_data)
#> Error in nlsModel(formula, mf, start, wts): singular gradient matrix at initial parameter estimates

Any help is greatly appreciated.

Solution

I assume you are given a vector of observations of length len as the ones plotted in your example, and you wish to identify k jumps and k jump sizes. (Or maybe I misunderstood you; but you have not really said what you want to achieve.) Below I will sketch a solution using Local Search. I start with your example data:

x <- seq(0, 10, by = 0.01)
y <- staircase(x,
               c(1,2,2,5),
               c(2,5,2,1)) + rnorm(length(x), mean = 0, sd = 0.2)

A solution is a list of positions and sizes of the jumps. Note that I use vectors to store these data, as it will become cumbersome to define variables when you have 20 jumps, say.

An example (random) solution:

k <- 5   ## number of jumps
len <- length(x)

sol <- list(position = sample(len, size = k),
            size = runif(k))

## $position
## [1]  89 236 859 885 730
## 
## $size
## [1] 0.2377453 0.2108495 0.3404345 0.4626004 0.6944078

We need an objective function to compute the quality of the solution. I also define a simple helper function stairs, which is used by the objective function. The objective function abs_diff computes the average absolute difference between the fitted series (as defined by the solution) and y.

stairs <- function(len, position, size) {
    ans <- numeric(len)
    ans[position] <- size
    cumsum(ans)
}

abs_diff <- function(sol, y, stairs, ...) {
    yy <- stairs(length(y), sol$position, sol$size)
    sum(abs(y - yy))/length(y)
}

Now comes the key component for a Local Search: the neighbourhood function that is used to evolve the solution. The neighbourhood function takes a solution and changes it slightly. Here, it will either pick a position or a size and modify it slightly.

neighbour <- function(sol, len, ...) {
    p <- sol$position
    s <- sol$size

    if (runif(1) > 0.5) {
        ## either move one of the positions ...
        i <- sample.int(length(p),  size = 1)
        p[i] <- p[i] + sample(-25:25, size = 1)
        p[i] <- min(max(1, p[i]), len)        
    } else {
        ## ... or change a jump size
        i <- sample.int(length(s), size = 1)
        s[i] <- s[i] + runif(1, min = -s[i], max = 1)
    }

    list(position = p, size = s)
}

An example call: here the new solution has its first jump size changed.

## > sol
## $position
## [1]  89 236 859 885 730
## 
## $size
## [1] 0.2377453 0.2108495 0.3404345 0.4626004 0.6944078
## 
## > neighbour(sol, len)
## $position
## [1]  89 236 859 885 730
## 
## $size
## [1] 0.2127044 0.2108495 0.3404345 0.4626004 0.6944078

I remains to run the Local Search.

library("NMOF")
sol.ls <- LSopt(abs_diff,
                list(x0 = sol, nI = 50000, neighbour = neighbour),
                stairs = stairs,
                len = len,
                y = y)

We can plot the solution: the fitted line is shown in blue.

plot(x, y)
lines(x, stairs(len, sol.ls$xbest$position, sol.ls$xbest$size),
      col = "blue", type = "S")