26
$\begingroup$

I'm trying to fit a line+exponential curve to some data. As a start, I tried to do this on some artificial data. The function is: $$y=a+b\cdot r^{(x-m)}+c\cdot x$$ It is effectively an exponential curve with a linear section, as well as an additional horizontal shift parameter (m). However, when I use R's nls() function I get the dreaded "singular gradient matrix at initial parameter estimates" error, even if I use the same parameters that I used to generate the data in the first place.
I've tried the different algorithms, different starting values and tried to use optim to minimise the residual sum of squares, all to no avail. I've read that a possible reason for this could be an over-parametrisation of the formula, but I don't think it is (is it?)
Does anyone have a suggestion for this problem? Or is this just an awkward model?

A short example:

#parameters used to generate the data
reala=-3
realb=5
realc=0.5
realr=0.7
realm=1
x=1:11 #x values - I have 11 timepoint data
#linear+exponential function
y=reala + realb*realr^(x-realm) + realc*x
#add a bit of noise to avoid zero-residual data
jitter_y = jitter(y,amount=0.2)
testdat=data.frame(x,jitter_y)

#try the regression with similar starting values to the the real parameters
linexp=nls(jitter_y~a+b*r^(x-m)+c*x, data=testdat, start=list(a=-3, b=5, c=0.5, r=0.7, m=1), trace=T)

Thanks!

$\endgroup$
2
  • 3
    $\begingroup$ Hint: look at the coefficient of $r^x$ (for a fixed $r$) and note that $b r^{-m} = \text{constant}$ has a one-dimensional family of solutions $(b,m)$ with $b = r^m \cdot \text{constant}$. $\endgroup$
    – whuber
    Commented Jul 14, 2011 at 17:35
  • 2
    $\begingroup$ This isn't an identified model, unless $b$ or $r$ are somehow constrained. I think requiring the $r \in (0,1)$ would do the job. $\endgroup$
    – Macro
    Commented Jul 14, 2011 at 19:02

2 Answers 2

21
$\begingroup$

I've got bitten by this recently. My intentions were the same, make some artificial model and test it. The main reason is the one given by @whuber and @marco. Such model is not identified. To see that, remember that NLS minimizes the function:

$$\sum_{i=1}^n(y_i-a-br^{x_i-m}-cx_i)^2$$

Say it is minimized by the set of parameters $(a,b,m,r,c)$. It is not hard to see that the the set of parameters $(a,br^{-m},0,r,c)$ will give the same value of the function to be minimized. Hence the model is not identified, i.e. there is no unique solution.

It is also not hard to see why the gradient is singular. Denote

$$f(a,b,r,m,c,x)=a+br^{x-m}+cx$$

Then

$$\frac{\partial f}{\partial b}=r^{x-m}$$

$$\frac{\partial f}{\partial m}=-b\ln rr^{x-m}$$

and we get that for all $x$

$$b\ln r\frac{\partial f}{\partial b}+\frac{\partial f}{\partial m}=0.$$

Hence the matrix

\begin{align} \begin{pmatrix} \nabla f(x_1)\\\\ \vdots\\\\ \nabla f(x_n) \end{pmatrix} \end{align}

will not be of full rank and this is why nls will give the the singular gradient message.

I've spent over a week looking for bugs in my code elsewhere till I noticed that the main bug was in the model :)

$\endgroup$
8
  • 2
    $\begingroup$ This is ages old I know but just wondering, does this mean nls cannot be used on models that are not identifiable? For example, a neural network? $\endgroup$
    – Count Zero
    Commented Jul 20, 2015 at 19:47
  • 1
    $\begingroup$ fat chance, I know, but could you break this down for less calc-remembering folks? :). also, what is the solution to the OP's problem, then? Give up and go home? $\endgroup$ Commented May 4, 2018 at 19:37
  • 3
    $\begingroup$ The solution to OP problem is to use one parameter instead of two, i.e. instead of $b\cdot r^{x-m}$ use $\beta \cdot r^x$. The $m$ parameter is superflous, since it merges into $\beta$, i.e. $\beta = b\cdot r^{-m}$. $\endgroup$
    – mpiktas
    Commented May 8, 2018 at 6:44
  • $\begingroup$ @CountZero, basically yes, usual optimisation methods would fail if the parameters are unindentified. Neural networks sidestep this problem however, by adding additional contraints and using other interesting tricks. $\endgroup$
    – mpiktas
    Commented May 8, 2018 at 6:47
  • $\begingroup$ I would think $\frac {\partial f}{\partial m} = -b \ln{r}\ r^{x-m}$ ? there is a lnr missing in the answer ? $\endgroup$
    – wiswit
    Commented Feb 19, 2019 at 5:54
18
$\begingroup$

The answers above are, of course, correct. For what its worth, in addition to the explanations given, if you are trying this on an artificial data set, according to the nls help page found at: http://stat.ethz.ch/R-manual/R-patched/library/stats/html/nls.html

R's nls wont be able to handle it. The help page specifically states:

Warning

Do not use nls on artificial "zero-residual" data.

The nls function uses a relative-offset convergence criterion that compares the numerical imprecision at the current parameter estimates to the residual sum-of-squares. This performs well on data of the form

y = f(x, θ) + eps

(with var(eps) > 0). It fails to indicate convergence on data of the form

y = f(x, θ)

because the criterion amounts to comparing two components of the round-off error. If you wish to test nls on artificial data please add a noise component, as shown in the example below.

So, no noise==no good for R's nls.

$\endgroup$
1
  • 2
    $\begingroup$ I am aware of this issue - hence using the "jitter" function to add some noise $\endgroup$
    – steiny
    Commented Aug 30, 2017 at 13:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.