question about nls

Actually, it likely won't matter where you start. The Gauss-Newton direction
is nearly always close to 90 degrees from the gradient, as seen by turning
trace=TRUE in the package nlmrt function nlxb(), which does a safeguarded
Marquardt calculation. This can be used in place of nls(), except you need
to put your data in a data frame. It finds a solution pretty
straightforwardly, though with quite a few iterations and function
evaluations.

Interesting observation but it does converge in 5 iterations with the
improved starting value whereas it fails due to a singular gradient
with the original starting value.

Lines <- "

+ x    y
+ 60 0.8
+ 80 6.5
+ 100 20.5
+ 120 45.9
+ "

DF <- read.table(text = Lines, header = TRUE)

# original starting value - singular gradient
nls(y ~ exp(a + b*x)+d,DF,start=list(a=0,b=0,d=1))

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

# better starting value - converges in 5 iterations
lm1 <- lm(log(y) ~ x, DF)
st <- setNames(c(coef(lm1), 0), c("a", "b", "d"))
nls(y ~ exp(a + b*x)+d, DF, start=st)

Nonlinear regression model
  model:  y ~ exp(a + b * x) + d
   data:  DF
      a       b       d
-0.1492  0.0342 -6.1966
 residual sum-of-squares: 0.5743

Number of iterations to convergence: 5
Achieved convergence tolerance: 6.458e-07

As Gabor indicates, using a start based on a good approximation is 
usually helpful, and nls() will generally find solutions to problems 
where there are such starts, hence the SelfStart methods. The Marquardt 
approaches are more of a pit-bull approach to the original 
specification. They grind away at the problem without much finesse, but 
generally get there eventually. If one is solving lots of problems of a 
similar type, good starts are the way to go. One-off (or being lazy), I 
like Marquardt.

It would be interesting to know what proportion of random starting 
points in some reasonable bounding box get the "singular gradient" 
message or other early termination with nls() vs. a Marquardt approach, 
especially as this is a tiny problem. This is just one example of the 
issue R developers face in balancing performance and robustness. The GN 
method in nls() is almost always a good deal more efficient than 
Marquardt approaches when it works, but suffers from a fairly high 
failure rate.


JN

On Fri, Mar 15, 2013 at 9:45 AM, Prof J C Nash (U30A) <nashjc at uottawa.ca> wrote:

Actually, it likely won't matter where you start. The Gauss-Newton direction
is nearly always close to 90 degrees from the gradient, as seen by turning
trace=TRUE in the package nlmrt function nlxb(), which does a safeguarded
Marquardt calculation. This can be used in place of nls(), except you need
to put your data in a data frame. It finds a solution pretty
straightforwardly, though with quite a few iterations and function
evaluations.

Interesting observation but it does converge in 5 iterations with the
improved starting value whereas it fails due to a singular gradient
with the original starting value.

Lines <- "

+ x    y
+ 60 0.8
+ 80 6.5
+ 100 20.5
+ 120 45.9
+ "

DF <- read.table(text = Lines, header = TRUE)

# original starting value - singular gradient
nls(y ~ exp(a + b*x)+d,DF,start=list(a=0,b=0,d=1))

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

# better starting value - converges in 5 iterations
lm1 <- lm(log(y) ~ x, DF)
st <- setNames(c(coef(lm1), 0), c("a", "b", "d"))
nls(y ~ exp(a + b*x)+d, DF, start=st)

Nonlinear regression model
   model:  y ~ exp(a + b * x) + d
    data:  DF
       a       b       d
-0.1492  0.0342 -6.1966
  residual sum-of-squares: 0.5743

Number of iterations to convergence: 5
Achieved convergence tolerance: 6.458e-07

I decided to follow up my own suggestion and look at the robustness of 
nls vs. nlxb. NOTE: this problem is NOT one that nls() would usually be 
applied to. The script below is very crude, but does illustrate that 
nls() is unable to find a solution in >70% of tries where nlxb (a 
Marquardt approach) succeeds.

I make no claim for elegance of this code -- very quick and dirty.

JN



debug<-FALSE
library(nlmrt)
x<-c(60,80,100,120)
y<-c(0.8,6.5,20.5,45.9)
mydata<-data.frame(x,y)
mydata
xmin<-c(0,0,0)
xmax<-c(8,8,8)
set.seed(123456)
nrep <- as.numeric(readline("Number of reps:"))
pnames<-c("a", "b", "d")
npar<-length(pnames)
# set up structure to record results
#  need start, failnls, parnls, ssnls, failnlxb, parnlxb, ssnlxb
tmp<-matrix(NA, nrow=nrep, ncol=3*npar+4)
outcome<-as.data.frame(tmp)
rm(tmp)
colnames(outcome)<-c(paste("st-",pnames[[1]],''),
	paste("st-",pnames[[2]],''),
	paste("st-",pnames[[3]],''),
         "failnls", paste("nls-",pnames[[1]],''),
	paste("nls",pnames[[1]],''),
	paste("nls-",pnames[[1]],''), "ssnls",
         "failnlxb", paste("nlxb-",pnames[[1]],''),
	paste("nlxb",pnames[[1]],''),
	paste("nlxb-",pnames[[1]],''), "ssnlxb")


for (i in 1:nrep){

cat("Try ",i,":\n")

st<-runif(3)
names(st)<-pnames
if (debug) print(st)
rnls<-try(nls(y ~ exp(a + b*x)+d,start=st, data=mydata), silent=TRUE)
if (class(rnls) == "try-error") {
    failnls<-TRUE
    parnls<-rep(NA,length(st))
    ssnls<-NA
} else {
    failnls<-FALSE
    ssnls<-deviance(rnls)
    parnls<-coef(rnls)
}
names(parnls)<-pnames
if (debug) {
   cat("nls():")
   print(rnls)
}
rnlxb<-try(nlxb(y ~ exp(a + b*x)+d,start=st, data=mydata), silent=TRUE)
if (class(rnlxb) == "try-error") {
    failnxlb<-TRUE
    parnlxb<-rep(NA,length(st))
    ssnlxb<-NA
} else {
    failnlxb<-FALSE
    ssnlxb<-rnlxb$ssquares
    parnlxb<-rnlxb$coeffs
}
names(parnls)<-pnames
if (debug) {
   cat("nlxb():")
   print(rnlxb)
   tmp<-readline()
   cat("\n")
   }
  solrow<-c(st, failnls=failnls, parnls, ssnls=ssnls,
      failnlxb=failnlxb, parnlxb, ssnlxb=ssnlxb)
   outcome[i,]<-solrow
} # end loop

cat("Proportion of nls  runs that failed = ",sum(outcome$failnls)/nrep,"\n")
cat("Proportion of nlxb runs that failed = 
",sum(outcome$failnlxb)/nrep,"\n")