1
$\begingroup$

I have a very large data set (roughly 11,000 points) that follow a roughly exponential curve with random variation. Here's a small sample of my data set:

Semi-exponential Graph

The underlying exponential function follows a curve a*b^x, not a*e^(b*x). In the case of the graph above, the equation is something like 0.05 * 1.195^x (I don't know the exact values)

My aim is to find the two parameters (a and b) that best fit the data. If the equation were of the form a*e^(b*x) then I could use standard linear regression techniques, but that's not the case.

So far I've taken the approach of computing the sum of squared errors (for all 11,000 data points) and attempting to minimize this error function. I've had varying degrees of success with this:

  • Using Excel I can set up two "input" cells to be my parameters a and b, add a column for "computed value" beside my data, add a column for "squared error", then add an "output" cell that sums this squared error column. I then use Solver to minimize the output cell by changing the input cells and it performs very well. On my raw data set I get a sum squared error around 48,000
  • Scipy's curve_fit utterly failed when passed my data set, giving a nonsensical answer. Scipy's minimize equally failed when passed an error function that computed the sum of squared errors but gave a meaningful message that overflow was encountered and so the desired error was not achieved due to precision loss. The error was in the range of 1e147
  • OpenOffice also has a Solver, like Excel, but it utterly failed and generated nonsense answers. The error was in the range of 1e9
  • Using both the genetic optimization and Powell optimization from optimization-js returned nonsense answers. The error was in the range of 1e9. I wasn't sure if I could take advantage of the L-BFGS or gradient descent options since I don't know how to compute the derivative for my error function
  • I tried to use liboptim but couldn't figure out how to get all of the dependencies working (Armadillo, LAPACK, etc)
  • I wrote my own very naive nonlinear solver which uses a pseudo-binary search for the first parameter and steps through all possible digits for the second parameter, stopping when it hits an inflection point. This seemed to work pretty well. It returned an error of 55,000 (not the 48,000 Excel got, but far better than I could do by hand)

In trying to research the mathematics involved in non-linear optimization so I can improve my naive optimizer, I keep stumbling upon recurring terms that I don't fully understand. Like whether a function is "Lipschitz", or whether the function is "convex".

My first question is: Given the definition of my error function (the sum of squared differences between an exponential function and a data set that is roughly exponential), what properties would my function have? Is it convex? Is it "Lipschitz"?

My second question is: Am I overdoing this? Is there an easier solution?

$\endgroup$
1
$\begingroup$

Your model is a pure exponential $$y=a\, b^x=a\, e^{x\log(b)}=a\, e^{cx}$$ but it is nonlinear with respect to its parameters; so you need some reasonable guesses to start.

Keeping your formulation, in a first step, linearize the model $$y=a\, b^x \implies \log(y)=\log(a)+x \log(b)=\alpha + \beta x$$ A first linear regression gives $\alpha$ and $\beta$ and then $a=e^{\alpha}$ and $b=e^{\beta}$. Now, start the nonlinear regression.

Edit

You could even reduce the problem to one equation in $b$ $$a=\frac{\sum_{i=1}^n y_i b^{x_i} } {\sum_{i=1}^n b^{2x_i} }$$ and then $$f(b)=\frac{\sum_{i=1}^n y_i b^{x_i} } {\sum_{i=1}^n b^{2x_i} }-\frac{\sum_{i=1}^n x_iy_i b^{x_i} } {\sum_{i=1}^n x_ib^{2x_i} }=0$$ SInce you have the estimate from $\beta$, even plotting will give you the result

$\endgroup$
6
  • $\begingroup$ Wait.... e^(x*log(b)) == b^x? You're blowing my mind right now. Either my college did me a massive disservice, or I've managed to forget something incredibly important after 12 years outside of a calculus classroom $\endgroup$
    – stevendesu
    Oct 4 '20 at 5:05
  • $\begingroup$ So since you pointed out that this COULD simplify to a linear equation, I attempted a linear regression (using regression-js and it came up with values of a=0.0463 and b=1.2054, which led to a sum of squared errors of 59122 for my data set -- very good, but not quite the 48000 I got from Excel (not even the 55000 I got with my naive nonlinear solver). Do you know why the squared error term wouldn't have been minimized with this method? Could it be an implementation issue in regression-js? $\endgroup$
    – stevendesu
    Oct 4 '20 at 5:22
  • 1
    $\begingroup$ @stevendesu. You must continue with the nonlinear regression since what is measured is $y$ and not $\log(y)$. $\endgroup$ Oct 4 '20 at 5:26
  • $\begingroup$ Using a and b from the linear regression as initial parameters to my nonlinear solver sped things up noticeable. Gradient descent completed in 6087ms instead of the previous ~15000ms (nearly a 3x gain). On a related note: I managed to get a gradient descent optimizer working (it required a gradient function which I figured out how to derive) $\endgroup$
    – stevendesu
    Oct 4 '20 at 5:33
  • $\begingroup$ @stevendesu. You do not need optimization if you use what I wrote in the edit. A root finding algorithm (Newton) with the $b$ from linear regression and numerical derivatives of $f(b)$ would suffice. $\endgroup$ Oct 4 '20 at 5:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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