Hi, All:
What tools exist for diagnosing singular gradient problems with
'nls'? Consider the following toy example:
DF1 <- data.frame(y=1:9, one=rep(1,9))
nlsToyProblem <- nls(y~(a+2*b)*one, DF1, start=list(a=1, b=1),
control=nls.control(warnOnly=TRUE))
Error in nlsModel(formula, mf, start, wts) :
singular gradient matrix at initial parameter estimates
This example is obviously stupid, but other singular gradient
problems are not so obvious.
If we transfer this problem to 'optim', we can get diagnostics
from an eigen analysis of the hessian:
dumfun <- function(x, y, one){
d <- y-(x[1]+2*x[2])*one
sum(d^2)
}
optimToyProblem <- optim(c(a=1, b=1), dumfun, hessian=TRUE,
y=DF1$y, one=DF1$one)
eigen(optimToyProblem$hessian, symmetric=TRUE)
$values
[1] 9.000000e+01 -7.105427e-10
$vectors
[,1] [,2]
[1,] 0.4472136 -0.8944272
[2,] 0.8944272 0.4472136
The smallest eigenvalue is essentially numerically zero relative
to the largest,confirming the 'singular gradient' message. The
corresponding eigenvector helps diagnose the problem: Adding (-0.9,
0.45)*z to any solution gives another equally good solution, for any z.
I've used this technique to diagnose many subtle convergence
problems with 'optim'. Are tools of this nature available for
'nls'?
Thanks,
Spencer Graves
Dear Spencer, I just saw your post. If the singular gradient happens during or after iteration one (that is, not at the initial estimates), then calling summary on the nls output would give standard error estimates on the parameters useful for diagnostics. You could also call chol2inv(xx$m$Rmat()) where xx is the object returned by nls to get an estimate of the inverse of the hessian; you could use this estimate to proceed with the diagnostics you were discussing. It would be possible (and in my opinion, desirable) to modify nls to return an object that contains the information above even if the singular gradient is at the intial estimates, but I don't think this can be accomplished without quite a few changes. You could also use nls.lm from the package minpack.lm to get a hessian estimate. By the way, I liked your idea from a while back (https://stat.ethz.ch/pipermail/r-devel/2008-April/048984.html) to do some diagnostics to identify problems with overparameterization automatically. best, Kate On Fri, 23 May 2008, Spencer Graves wrote:> Hi, All: > > What tools exist for diagnosing singular gradient problems with > 'nls'? Consider the following toy example: > > DF1 <- data.frame(y=1:9, one=rep(1,9)) > nlsToyProblem <- nls(y~(a+2*b)*one, DF1, start=list(a=1, b=1), > control=nls.control(warnOnly=TRUE)) > Error in nlsModel(formula, mf, start, wts) : > singular gradient matrix at initial parameter estimates > > This example is obviously stupid, but other singular gradient > problems are not so obvious. > > If we transfer this problem to 'optim', we can get diagnostics > from an eigen analysis of the hessian: > > dumfun <- function(x, y, one){ > d <- y-(x[1]+2*x[2])*one > sum(d^2) > } > optimToyProblem <- optim(c(a=1, b=1), dumfun, hessian=TRUE, > y=DF1$y, one=DF1$one) > eigen(optimToyProblem$hessian, symmetric=TRUE) > $values > [1] 9.000000e+01 -7.105427e-10 > > $vectors > [,1] [,2] > [1,] 0.4472136 -0.8944272 > [2,] 0.8944272 0.4472136 > > The smallest eigenvalue is essentially numerically zero relative > to the largest,confirming the 'singular gradient' message. The > corresponding eigenvector helps diagnose the problem: Adding (-0.9, > 0.45)*z to any solution gives another equally good solution, for any z. > > I've used this technique to diagnose many subtle convergence > problems with 'optim'. Are tools of this nature available for 'nls'? > > Thanks, > Spencer Graves > > ______________________________________________ > R-help at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. >
On Sun, May 25, 2008 at 9:32 PM, Katharine Mullen <kate at few.vu.nl> wrote:> Dear Spencer, > > I just saw your post. > > If the singular gradient happens during or after iteration one (that is, > not at the initial estimates), then calling summary on the nls output > would give standard error estimates on the parameters useful for > diagnostics. You could also call chol2inv(xx$m$Rmat()) where xx is the > object returned by nls to get an estimate of the inverse of the hessian; > you could use this estimate to proceed with the diagnostics you were > discussing.Try this:> library(nls2) > DF1 <- data.frame(y=1:9, one=rep(1,9)) > xx <- nls2(y~(a+2*b)*one, DF1, start = c(a=1, b=1), algorithm = "brute-force") > eigen(chol2inv(xx$m$Rmat()))$values [1] 5.070602e+31 0.000000e+00 $vectors [,1] [,2] [1,] -0.8944272 -0.4472136 [2,] 0.4472136 -0.8944272
Hi, Gabor: Thanks very much. Spencer Gabor Grothendieck wrote:> On Sun, May 25, 2008 at 9:32 PM, Katharine Mullen <kate at few.vu.nl> wrote: > >> Dear Spencer, >> >> I just saw your post. >> >> If the singular gradient happens during or after iteration one (that is, >> not at the initial estimates), then calling summary on the nls output >> would give standard error estimates on the parameters useful for >> diagnostics. You could also call chol2inv(xx$m$Rmat()) where xx is the >> object returned by nls to get an estimate of the inverse of the hessian; >> you could use this estimate to proceed with the diagnostics you were >> discussing. >> > > Try this: > > >> library(nls2) >> DF1 <- data.frame(y=1:9, one=rep(1,9)) >> xx <- nls2(y~(a+2*b)*one, DF1, start = c(a=1, b=1), algorithm = "brute-force") >> eigen(chol2inv(xx$m$Rmat())) >> > $values > [1] 5.070602e+31 0.000000e+00 > > $vectors > [,1] [,2] > [1,] -0.8944272 -0.4472136 > [2,] 0.4472136 -0.8944272 > > ______________________________________________ > R-help at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. >