I have a table of data with input values and target values, and I was tasked to do something quite peculiar. I was tasked to run a sort of exponential regression and report back with the coefficients and the series of exponents. Let me explain:

Vanilla Linear Regression will minimize sum of squared error to create a linear combination in the format of $[a_1x_1 + a_2x_2...]$ where $a_i$ represents the series of coefficients and $x_i$ represents the different input features. This means that solving for the coefficients on linear regression is fairly easy given enough linearly independent data points.

I was given a dataset, and I was tasked to create an equation in the form of $[a_1x_1^{n_1} + a_2x_2^{n_2}...]$ where $a_i$ is the coefficients, $x_i$ is the input features, and $n_i$ signifies an exponential value.

I believe that since $a_i$ and $n_i$ need to be solved for, there are too many variables involved to go about this problem in a simple manner.

Is something like this even something like this even solvable? If so, how would you go about it statistically, and how would you go about it in python?

p.s. The task also says that all exponents, $n_i$, need to be between 0 and 1, so that might make it a bit easier.

p.s.s. I know that $[a_1x_1^{n_1} + a_2x_2^{n_2}...]$ is not linear unless all $n_i$ variables are either 0 or 1 but I wrote "linear regression" in the title because the two general formulas are similar.

  • $\begingroup$ The condition to be linear regression is $n_i$ be a known constant. $\endgroup$ – user158565 Aug 8 '19 at 17:46
  • $\begingroup$ If we approach this sequentially, you could find the n exponents that makes the relationships for each x with y most linear via a Box-Cox transformation, and the run a linear regression using the transformed x variables. But I'm confused as to why the n must be between 0 and 1. You might want to wait and see if you get a better answer. $\endgroup$ – zbicyclist Aug 8 '19 at 18:22
  • $\begingroup$ This is called fractional polynomial regression. It can be fit using exploratory methods as suggested by @zbicyclist or, more formally, by a "backward elimination" method that focuses on one variable at a time. If you truly need all exponents between $0$ and $1$ that's an unwelcome complication: it's simpler to find out the best values of the exponents so you don't have to deal with any constraints. $\endgroup$ – whuber Aug 8 '19 at 21:16

Here is example code that I think should what you describe, using the equation "Z = (a1 * X1^n1) + (a2 * X2^n2)". This code uses the "brick wall" technique to keep the parameters n1 and n2 within bounds of 0.0 and 1 as noted in the comments in func().

fitted prameters:

a1 = 1.15465275488

n1 = 0.947251333112

a2 = 0.00698288550135

n2 = 0.999999491568

RMSE: 0.189049737197

R-squared: 0.99528439111

import numpy, scipy, scipy.optimize

def func(data, a1, n1, a2, n2):
    # if n1 or n2 is less than zero or greater than
    # one, return a large error - this will act as
    # a "brick wall" that the curve fit cannot overcome,
    # and this will act to keep the parameters within these
    # bounds. Note that the initial parameters estimates for
    # n1 and n2 (set below) must be within these bounds.
    if n1 < 0.0 or n1 > 1.0:
        return 1.0E300
    if n2 < 0.0 or n2 > 1.0:
        return 1.0E300

    x1 = data[0]
    x2 = data[1]
    return (a1 * x1**n1) + (a2 * x2**n2)

if __name__ == "__main__":
    x1Data = numpy.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0])
    x2Data = numpy.array([11.0, 12.1, 13.0, 14.1, 15.0, 16.1, 17.0, 18.1, 90.0])
    zData = numpy.array([1.1, 2.2, 3.3, 4.4, 5.5, 6.6, 7.7, 8.0, 9.9])

    data = [x1Data, x2Data, zData]

    # these are a1, n1, a2, n2 per the function definition
    # n1 and n2 are set within the bounds of 0.0 and 1.0
    initialParameters = [1.0, 0.5, 1.0, 0.5]

    # here a non-linear surface fit is made with scipy's curve_fit()
    fittedParameters, pcov = scipy.optimize.curve_fit(func, [x1Data, x2Data], zData, p0 = initialParameters)

    print('fitted prameters:')
    print('a1 =', fittedParameters[0])
    print('n1 =', fittedParameters[1])
    print('a2 =', fittedParameters[2])
    print('n2 =', fittedParameters[3])

    modelPredictions = func(data, *fittedParameters) 

    absError = modelPredictions - zData

    SE = numpy.square(absError) # squared errors
    MSE = numpy.mean(SE) # mean squared errors
    RMSE = numpy.sqrt(MSE) # Root Mean Squared Error, RMSE
    Rsquared = 1.0 - (numpy.var(absError) / numpy.var(zData))
    print('RMSE:', RMSE)
    print('R-squared:', Rsquared)

Fitting to experimental data the function : $$y(x)=b\:x^{p}+c\:x^{q}$$ for the approximate values of the 4 parameters $p,q,b,c$ is easy thanks to the method described pp.71-72 (with numerical example p.73) in the paper https://fr.scribd.com/doc/14674814/Regressions-et-equations-integrales

The case of the 5 parameters function : $$y(x)=a+b\:x^{p}+c\:x^{q}$$ is shown in https://math.stackexchange.com/questions/3306953/sum-of-exponential-growth-and-decay/3317563#3317563

The extension to more parameters is theoretically possible on the same way. But too many parameters increase the deviations in the numerical calculus, making the result less and less accurate, up to be not acceptable in practice. This is discussed p.74 in the referenced paper.

  • $\begingroup$ Your links relate to a sum of exponential functions of $x$, not to a sum of powers of $x$. $\endgroup$ – Sextus Empiricus Nov 24 '20 at 9:30
  • $\begingroup$ The method which you describe on pp. 71-72 for fitting a sum of exponentials seems to be the integration method by Tittelbach and Helmrich also described on Cross Validated here: stats.stackexchange.com/a/428847 $\endgroup$ – Sextus Empiricus Nov 24 '20 at 9:35
  • $\begingroup$ @Sextus Empiricus. The case of sum of powers of $x$ is treated in the paper referenced in my answer. The method is not iterative and doesn't requires initial values, that is no need for guess or preliminary calculus. $\endgroup$ – JJacquelin Nov 24 '20 at 10:12
  • $\begingroup$ You mention two sources but neither of them provide a (direct) solution to the problem. The pages 71 and 72 of your reference and the linked maths.stackexchange speak about a sum of exponential terms $e^{a_i x}$ not a sum of power terms $x^{a_i}$. Maybe your reference also speaks about a sum of powers but in that case you might need to refer to different page numbers. $\endgroup$ – Sextus Empiricus Nov 24 '20 at 10:31
  • $\begingroup$ I see now that I missed the case of powers on page 72. This is what makes link-based answers and referring to other sources without providing a small context difficult for this Q&A format. It is especially the link to the math.stackexchange that is confusing because it does not mention powers at all. $\endgroup$ – Sextus Empiricus Nov 24 '20 at 10:48

If $x_i$ are the same

If you transform the $x$ values like:

$$t = \log(x)$$

Then your formula changes into

$$y = a_1 \exp(n_1t) + a_2 \exp(n_2t) + \dots$$

This you can solve with one of many methods for sums of exponential functions.

If $x_i$ are not the same

The above methods are for the case when the features $x_i$ are the same feature. E.g. a time variable. Your case might be a bit more difficult.

Solving $$y = a_1 \exp(n_1t) + a_2 \exp(n_2t) + \dots$$ is easier than $$y = a_1 \exp(n_1t_1) + a_2 \exp(n_2t_2) + \dots$$

For the latter, there is no easy way to find a good starting value. However, finding the solution can be made slightly easier by the method described below.

Segregated solving

You can solve it in a segregated fashion by considering that the function is a separable non-linear function and can be divided in a linear part and a non-linear part. When we try to minimize the loss function $\mathcal{L}(a_1,a_2,t_1,t_2)$ we can eliminate the coefficients $a_1$ and $a_2$ by considering them as the solution for an OLS fit with $t_1$ and $t_2$ fixed. Then these coefficients $a_1$ and $a_2$ are simple linear functions of $t_1$ and $t_2$ and we fill them into the equation for the loss function $\mathcal{L}(f_1(t_1,t_2),f_2(t_1,t_2),t_1,t_2)$ which is now effectively a function of only two parameters $\mathcal{L}(t_1,t_2)$. In this question an example is given for this method: How to compute the gradient for a seperable nonlinear least squares problem?

The segregated solving will optimize the algorithm to find the minimum, but you still need to find a starting point for $t_1$ and $t_2$. This is a difficult problem as the cost function is not convex. The example below (from the previously linked question) shows how the solution may diverge away from the global minimum.

non convex example


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