# Exponential fitting - isn't actually a BUG there?

There are many questions on this site about wrong exponential fitting in Mathematica but no one considers this well-known problem as a potential bug. Usually people suggest well-known workarounds allowing to obtain expected results but they don't investigate the source of the problem. My recent attempt to understand the situation lead me to conclusion that there is a bug in FindFit and FindMinimum related to exponential fitting.

Let us take a simple dataset and try to fit it to the exponential model:

data = {{0, 20}, {20, 10}, {40, 5}, {60, 2.5}, {80, 1.25}};
fit1 = FindFit[data, a*Exp[b*x], {a, b}, x, EvaluationMonitor :> Print[{a, b}]]


{1.,1.}

{0.,1.}

 {a -> 0., b -> 1.}


Isn't the observed behavior already strange? Only two points are taken and then the last of them is returned as the minimum without any messages! It is easy to see that the returned solution is completely wrong and the correct solution can be obtained with Method -> NMinimize:

fit2 = FindFit[data, a*Exp[b*x], {a, b}, x, Method -> NMinimize]
Show[ListPlot[data], Plot[Evaluate[a*Exp[b*x] /. {fit1, fit2}], {x, 0, 80}]]

{a -> 20., b -> -0.0346574}


But WHY the default method (which is "LevenbergMarquardt") gives wrong result? Maybe there actually is a local minimum with zero gradient? Let us check the gradient of the Sum of Squared Residuals (which is the function minimized by FindFit; SSR below) at the points checked by FindFit:

Clear[a, b, x]
model[x_] := a*Exp[b*x];
resudual[{x_, y_}] := y - model[x];
residuals = Table[resudual[point], {point, data}];
SSR = Total[residuals^2];
gradientVector /. {{a -> 1., b -> 1.}, {a -> 0., b -> 1.}}

{{6.1397*10^69, 4.91176*10^71}, {-1.38516*10^35, 0.}}


The gradient is HUGE in absolute value at both points checked by FindFit but it returns the last of them as the solution without any messages! FindMinimum behaves identically:

FindMinimum[SSR, {a, b}]

{532.813, {a -> 0., b -> 1.}}


Can anyone explain this behavior (I'm not looking for a workaround, I'm looking for an explanation)? Is this a bug?

• Good question...(and an upvote). – Daniel Lichtblau Jan 27 '16 at 19:50
• Not an explanation you're asking for, but just an observation: Giving b a starting value of -1 yields the correct answer. – Sjoerd C. de Vries Jan 27 '16 at 21:31
• @DanielLichtblau Please see my answer. Having access to the sources, now it should be easy to figure out what happens. – Alexey Popkov Jan 28 '16 at 10:23
• It might help someone. Not my particular area. – Daniel Lichtblau Jan 28 '16 at 15:48

This looks like a numerical precision issue. Various approaches that precisely address this, all yield the same, correct solution:

Scaling of the x values:

data = {{0, 20}, {20, 10}, {40, 5}, {60, 2.5}, {80, 1.25}};
{xmin, xmax} = MinMax[data[[All,1]]];
data2 = {Rescale[data[[All, 1]]], data[[All, 2]]} // Transpose;
fit1 = FindFit[data2, a*Exp[b*x], {a, b}, x];
Show[ListPlot[data],
Plot[a*Exp[b*Rescale[x, {xmin, xmax}]] /. fit1, {x, 0, 80}]]


Exact numbers:

data = {{0, 20}, {20, 10}, {40, 5}, {60, 25/10}, {80, 125/100}};
fit1 = FindFit[data, a*Exp[b*x], {a, {b, -1}}, x]
Show[ListPlot[data], Plot[Evaluate[a*Exp[b*x] /. fit1], {x, 0, 80}]]
(* {a -> 20., b -> -0.0346574} *)


Arbitrary precision numbers:

data = {{0, 20}, {20, 10}, {40, 5}, {60, 2.5500}, {80, 1.25500}};
fit1 = FindFit[data, a*Exp[b*x], {a, {b, -1}}, x]

(*{a -> 19.9999999999999999999999999999999999999999999999999999999999999\
9999999999999999999999999999999999999999999999999999999999999999999999\
9999999999999999999999999999999999999999999999999999999999999999999999\
9999999999999999999999999999999999999999999999999999999999999999999999\
9999999999999999999999999778047394451665874718796271036850250884181562\
4585573342662728816721398861154072988885262039061332356855851201362915\
7164434036094129659039810884391575813515277252575297641755499155092912\
05514394088532664,
b -> -0.0346573590279972654708616060729088284037750067180127627060340\
0047466968109848473578029316634982093437710007405102853428668427601178\
7906527851633537581753798096536378541418571759515351931194583673556167\
5057682248977619560237586340787466032577762367069762941475226503547663\
1833213270521195789074760212376423504358288082923157416807773455937594\
8730959223814463379866456988534715801639211694432455726463474594405763\
7863828006795995086471222261587769767698218890439528263697455009896183\
99973828035023378735} *)

• It's sufficient in the last case just to add WorkingPrecision -> $MachinePrecision. It is b that needs to have the arbitrary precision, rather than the inputs, it seems. For machine-precision b it will fail even with increased WorkingPrecision. – Oleksandr R. Jan 28 '16 at 3:02 • Please see my answer. – Alexey Popkov Jan 28 '16 at 10:02 • The "Exact numbers" example doesn't work for me unless a helpful starting value is given for b. Also on my computer, giving almost any value for WorkingPrecision works (even lower than$MachinePrecision). But "WorkingPrecision->MachinePrecision" does not work. – Simplex Feb 3 '16 at 7:53
• @simplex are you sure that you have converted all numbers to exact numbers? I find it surprising that there should be a difference between our computers in this regard as the arbitrary precision routines are supposedly machine independent. – Sjoerd C. de Vries Feb 3 '16 at 7:56
• @SjoerdC.deVries, if I copy your code it works, but if I edit it fit1 = FindFit[data, a*Exp[b*x], {a, b}, x] then it does not. Does the edited version work for you? – Simplex Feb 3 '16 at 8:10

I am one of the many that does not consider this a bug as it is shared by probably the majority of statistical packages. (But there certainly is room for improvement in terms of warnings that could be given.)

Iterative procedures work great if the starting values are close enough to the final values. When one doesn't know the final result this requirement can be a bit daunting.

However, there are things that one can do to reduce the number of such failures.

(1) Find out what starting values are used by the fitting function. For this function Mathematica uses 1.0 as the starting value for all parameters. Consider the part of the function Exp[b*x]. With b = 1.0 and x = 80. Exp[1.0*80] = 5.54062*10^34. At the other end of the data there is x = 0 which gives Exp[1.0*0] = 1.0. This range of numerical values is bound to cause trouble.

(2) Standardize the independent variables. You can (a) subtract the mean and divide by the standard deviation or range, or (b) subtract the minimum value from each value if there are numerous superflous digits (as in geographical UTM units), or (c) divide by a power of 10 to get all values between 0 and 10. (And one should determine that doing so doesn't change the basic equation.)

One should also consider whether linearizing the equation might help in numerical stability and have the distribution of the residuals match the assumptions. (In other words, fitting y = logA + b * x + error.)

For this particular dataset one could divide the x values by 10 and everything works:

fit2 = FindFit[data, a*Exp[b*x/10], {a, b}, x,
EvaluationMonitor :> Print[{a, b}]]
(* {1.,1.}
{0.041318727177279246,0.9948817357756237}
{0.04283227733515496,0.8664735714509461}
{0.08908837807900087,0.6097184994465492}
{0.3761606155865537,0.09730611472542439}
{5.268929789483257,-0.796208595732955}
{13.173242486336434,-0.5173835399541626}
{19.927460536660046,-0.17029501801604285}
{19.503033453296695,-0.2801530807616852}
{19.9570169999364,-0.33847226843982786}
{19.999529625704145,-0.3464506437357793}
{19.999999896325985,-0.3465735609902636}
{19.999999999999993,-0.3465735902799709} *)


Update

Maybe to show what is going on (which I'm still claiming is a scaling issue) is the following:

Manipulate[
(* Determine fit *)
fit = Reap[
FindFit[data, a*Exp[b*x/k], {{a, aa}, {b, bb}}, x,
EvaluationMonitor :> Sow[{a, b}]]];

(* Adjust estimate of b to be on original scale *)
fit[[2, 1, All, 2]] = fit[[2, 1, All, 2]]/k;

(* Plot results *)
ListPlot[{{{20, -0.0346754}}, Prepend[fit[[2, 1]], {aa, bb}],
Prepend[fit[[2, 1]], {aa, bb}], {{a, b/k}} /. fit[[1]]},
PlotStyle -> {{Red, PointSize[0.05]}, {Blue,
PointSize[0.02]}, {Blue, PointSize[0.02]}, {Green,
PointSize[0.02]}},
Joined -> {False, True, False}, AxesOrigin -> {-1, 0},
PlotRange -> {{-1, 22}, {ymin, ymax}},
Frame -> True,
FrameLabel -> {Style["a", Large, Bold, Black],
Rotate[Style["b", Large, Bold, Black], -\[Pi]/2]}],

(* Controls *)
{{aa, 1, "a"}, -5, 22, Appearance -> "Labeled"},
{{bb, 1, "b"}, -0.5, 2, Appearance -> "Labeled"},
{{k, 2, "Divisor"}, 0.1, 100, Appearance -> "Labeled"},
{{ymax, 1.2}, 0, 1.5, Appearance -> "Labeled"},
{{ymin, -0.3}, -1, 0.1, Appearance -> "Labeled"},
TrackedSymbols :> {aa, bb, k, ymin, ymax},
Initialization :> (data = {{0, 20}, {20, 10}, {40, 5}, {60,
2.5}, {80, 1.25}})]


I believe this shows that without scaling of the x values, the search algorithm is dominated by looking for good values of b and essentially ignoring a until b is found. A more stable search occurs when successive points move a bit in both a and b directions.

Here is where the scaling factor (i.e., what we divide b by) is 2:

Now with the scaling factor being 100:

• Note that I've essentially just expanded on @SjoerdC.deVries section on "Scaling of the X values". – JimB Jan 28 '16 at 0:07
• It will be correct if you put 0.6 as an initial guess for b (Exp[80 b] == 7.0*^20), not converge for 0.7 (2.0*^24), and not even start to optimize for 0.8 (6.2*^27). The list of sampling points is interesting in these cases; it looks like it is having difficulty establishing a gradient with respect to a. – Oleksandr R. Jan 28 '16 at 3:00
• For whatever it's worth, PROC NLIN in $SAS$ gives 0 and 1 for the estimates of a and b, respectively, when using a*exp(b*x) but gives the correct results for a*exp(b*x/10). – JimB Jan 28 '16 at 4:00
• @OleksandrR. From the other side, the observed sequence of points checked by FindFit ({1., 1.}, then {0., 1.}) does not look like being produced by the Levenberg-Marquardt algorithm but looks like a built-in sequence to check before starting actual Levenberg-Marquardt. – Alexey Popkov Jan 28 '16 at 4:53
• @OleksandrR. When we specify starting value for b as 0.8 or more the algorithm checks 1. and 0. for a and it also seems to be a built-in sequence, not Levenberg-Marquardt. So the bug is somewhere in "preliminary" code which decides do not start Levenberg-Marquardt searching at all. – Alexey Popkov Jan 28 '16 at 5:02

In the comments I supposed that the reason why the algorithm stops is obtaining zero derivative with respect to b at the point {a, b} = {0., 1.} returned as the minimum. Let us check this statement by perturbing the actual Jacobian a little bit in order to make all the derivatives non-zero.

Proceeding from the code in the question, we find the Jacobian and feed it to FindMinimum as described here:

jacobianMatrix = D[residuals, {{a, b}}];

FindMinimum[, {a, b},
Method -> {"LevenbergMarquardt", "Residual" :> Sqrt[2] residuals,
"Jacobian" :> (Sqrt[2] jacobianMatrix)}]

{532.813, {a -> 0., b -> 1.}}


The result is exactly the same as with automatic Jacobian:

% === FindMinimum[SSR, {a, b}]

True


Now we add a tiny constant perturbation to the Jacobian:

FindMinimum[, {a, b},
Method -> {"LevenbergMarquardt", "Residual" :> Sqrt[2] residuals,
"Jacobian" :> (Sqrt[2] jacobianMatrix - 10^-150)}]


FindMinimum::cvmit: Failed to converge to the requested accuracy or precision within 100 iterations. >>

{531.25, {a -> 3.37441*10^-35, b -> 0.994968}}


The situation has moved forward! Let us increase the number of iterations:

FindMinimum[SSR, {a, b},
Method -> {"LevenbergMarquardt", "Residual" :> Sqrt[2] residuals,
"Jacobian" :> (Sqrt[2] jacobianMatrix - 10^-150)}, MaxIterations -> 10000]

{8.39587*10^-20, {a -> 20., b -> -0.0346574}}


We have got the correct solution!

Now let us investigate it a little deeper and create a "smarter" perturbation which will be non-zero only at the second evaluation of the Jacobian:

counter = 0;
myJacobian[count_] := If[count == 2, jacobianMatrix + 10^-150, jacobianMatrix];

FindMinimum[SSR, {a, b},
Method -> {"LevenbergMarquardt", "Residual" :> Sqrt[2] residuals,
"Jacobian" :> (Sqrt[2] myJacobian[++counter])}, MaxIterations -> 10000]

{5.29677*10^-25, {a -> 20., b -> -0.0346574}}


It is easy to check that if we make the perturbation non-zero at any other evaluation except the second, we obtain incorrect solution. So it is the second value of the Jacobian what is the stumbling block for the algorithm:

jacobianMatrix /. {a -> 0., b -> 1.}

{{-1, 0}, {-4.85165*10^8, 0.}, {-2.35385*10^17, 0.}, {-1.14201*10^26,
0.}, {-5.54062*10^34, 0.}}


One can see that all derivatives with respect to b are zero (as expected for a = 0).

Now let us make all the derivatives with respect to b zero at some arbitrary step of the algorithm (but at the second step we will keep the perturbation):

counter = 0;
myJacobian[count_] :=
Which[count == 2, jacobianMatrix + 10^-150, count == 10,
jacobianMatrix*ConstantArray[{1, 0}, 5], True, jacobianMatrix];

FindMinimum[SSR, {a, b},
Method -> {"LevenbergMarquardt", "Residual" :> Sqrt[2] residuals,
"Jacobian" :> (Sqrt[2] myJacobian[++counter])}, MaxIterations -> 10000]
counter

{531.25, {a -> 1.1088*10^-31, b -> 0.89375}}

11


Experimentation shows that the algorithm always stops after the step when it gets the Jacobian with zero derivatives with respect to b (but it won't stop if at least for one data point this derivative is non-zero). From the other side, if we make at some step zero derivatives with respect to a, it does not stop the algorithm:

counter = 0;
myJacobian[count_] :=
Which[count == 2, jacobianMatrix + 10^-150, count == 1000,
jacobianMatrix*ConstantArray[{0, 1}, 5], True, jacobianMatrix];

FindMinimum[SSR, {a, b},
Method -> {"LevenbergMarquardt", "Residual" :> Sqrt[2] residuals,
"Jacobian" :> (Sqrt[2] myJacobian[++counter])}, MaxIterations -> 10000]
counter

{1.97215*10^-31, {a -> 20., b -> -0.0346574}}

4119


Of course I still have not explained what is so special about zero derivatives with respect to b. I still feel this behavior as a bug or at least significant imperfection of current implementation.