Fit sum of exponentials

Question

I am looking for an algorithm that fits a sum of exponential. for example I have something like this: $$y(x)=ae^{−bx}+c+de^{-fx}+h$$ and I want to find a,b,c,d,f and h values. Of this sum I have only the x and y values that belong to the curve represented by the previous model. Is there any paper that explains that?

Answer 1

Edit: some people have asked me about the license of the code below. The license is MIT as specfified in this Github repository.

---

I would like to add another solution, after the useful discussion with @JJacquelin. I reworked my previous solution based on linear algebra and derivatives, to work with integrals instead. It also works for $n$ exponentials.

Assuming the data set $(y, x)$;

Compute the $n$ integrals of $y$ ; $\int y, \int^2{y}, \cdots , , \int^{n-1}{y} , \int^n{y}$.

Compute the $n-1$ powers of the independent variable $x$ ; $x, x^2, \cdots , x^{n-1}$.

Solve linear least squares problem for $[a_{1}, a_{2}, \cdots, a_{n-1}, a_{n}, k_{n-1}, k_{n-2} , \cdots , k_{2}, k_{1}, k_{0}]$ in:

$$ y = \left[ \int y , \int^2{y} , \cdots, \int^{n-1}{y}, \int^n{y} , x^{n-1}, x^{n-2}, \cdots , x^{2}, x, 1 \right] \left[ \begin{array}{c} a_{1} \\ a_{2} \\ \vdots \\ a_{n-1} \\ a_{n} \\ k_{n-1} \\ k_{n-2} \\ \vdots \\ k_{2} \\ k_{1} \\ k_{0} \end{array} \right] $$

$$ y = Y A \\ Y^T Y A = Y^T y \\ A = (Y^T Y)^{-1}Y^T y $$

Then use the first $n$ elements of $A$ to form the transfer matrix of the integral equation:

$$ \bar{A} = \left[ \begin{array}{cccccc} a_1 & a_2 & a_3 & \cdots & a_{n-1} & a_{n} \\ 1 & 0 & 0 & \cdots & 0 & 0 \\ 0 & 1 & 0 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & 1 & 0 \end{array} \right] $$

Compute the eigenvalues of $\bar{A}$ to obtain the solution of the characteristic polynomial.

The eigenvalues $\lambda _n$ are then the values of the exponentials in:

$$ y = p_1 e^{\lambda _1 x} + p_2 e^{\lambda _2 x} + \cdots + p_{n-1} e^{\lambda _{n-1} x} + p_n e^{\lambda _n x} $$

To obtain the missing $p_n$, calculate numerically the exponential terms $e^{\lambda _n}$ and solve the least squares problem for $[p_1, p_2, \cdots, p_{n-1}, p_n ]$ in

$$ y = [e^{\lambda _1 x}, e^{\lambda _2 x} , \cdots, e^{\lambda _{n-1} x}, e^{\lambda _n x}] \left[ \begin{array}{c} p_1 \\ p_2 \\ \vdots \\ p_{n-1} \\ p_n \end{array} \right] $$

$$ y = X P \\ X^T X P = X^T y \\ P = (X^T X)^{-1}X^T y $$

Now you have your $n$ exponential model parameters $p_n$ and $\lambda _n$.

This solution based on integrals, looks much more similar to @JJacquelin's solutions for cases $n = 2$ and $n = 3$, but now generalized for any $n$.

Here is the Matlab code for a test case with $n = 4$ (you can test it online here):

clear all;
clc;
% get data
dx = 0.02;
x  = (dx:dx:1.5)';
y  =  5*exp(0.5*x) + 4*exp(-3*x) + 2*exp(-2*x) - 3*exp(0.15*x);

% calculate integrals
iy1 = cumtrapz(x, y);
iy2 = cumtrapz(x, iy1);
iy3 = cumtrapz(x, iy2);
iy4 = cumtrapz(x, iy3);

% get exponentials lambdas
Y = [iy1, iy2, iy3, iy4, x.^3, x.^2, x, ones(size(x))];
A = pinv(Y)*y;

lambdas = eig([A(1), A(2), A(3), A(4); 1, 0, 0, 0; 0, 1, 0, 0; 0, 0, 1, 0]);
lambdas
%lambdas =
%  -2.9991
%  -1.9997
%   0.5000
%   0.1500

% get exponentials multipliers
X = [exp(lambdas(1)*x), exp(lambdas(2)*x), exp(lambdas(3)*x), exp(lambdas(4)*x)];
P = pinv(X)*y;
P
%P =
%   4.0042
%   1.9955
%   4.9998
%  -2.9996

Variations

Additive constant term

If it is desired to fit a sum of exponentials plus constant $k$:

$$ y = k + p_1 e^{\lambda _1 x} + p_2 e^{\lambda _2 x} + \cdots + p_{n-1} e^{\lambda _{n-1} x} + p_n e^{\lambda _n x} $$

Just add $x^n$ to the first least squares problem and a constant to the second least squares problem:

clear all;
clc;
% get data
dx = 0.02;
x  = (dx:dx:1.5)';
y  =  -1 + 5*exp(0.5*x) + 4*exp(-3*x) + 2*exp(-2*x);

% calculate integrals
iy1 = cumtrapz(x, y);
iy2 = cumtrapz(x, iy1);
iy3 = cumtrapz(x, iy2);

% get exponentials lambdas
Y = [iy1, iy2, iy3, x.^3, x.^2, x, ones(size(x))];
A = pinv(Y)*y;

lambdas = eig([A(1), A(2), A(3); 1, 0, 0; 0, 1, 0]);
lambdas
%lambdas =
%  -2.9991
%  -1.9997
%   0.5000

% get exponentials multipliers
X = [ones(size(x)), exp(lambdas(1)*x), exp(lambdas(2)*x), exp(lambdas(3)*x)];
P = pinv(X)*y;
P
%P =
%  -0.9996
%   4.0043
%   1.9955
%   4.9999

Sine and Cosine

If it is desired to fit a sum of exponentials plus a $cosine$ or a $sine$, the algorithm is exactly the same. Sines and cosines manifest themselves as a pair of complex conjugate eigen-values of $\hat{A}$.

For example, if we fit one exponential and one cosine:

clear all;
clc;
% get data
dx = 0.02;
x  = (dx:dx:1.5)';
y  =  5*exp(0.5*x) + 4*cos(0.15*x);

% calculate integrals
iy1 = cumtrapz(x, y);
iy2 = cumtrapz(x, iy1);
iy3 = cumtrapz(x, iy2);

% get exponentials lambdas
Y = [iy1, iy2, iy3, x.^2, x, ones(size(x))];
A = pinv(Y)*y;
%lambdas =
%   0.5000 +      0i
%  -0.0000 + 0.1500i
%  -0.0000 - 0.1500i

lambdas = eig([A(1), A(2), A(3); 1, 0, 0; 0, 1, 0]);
lambdas

% get exponentials multipliers
X = [exp(lambdas(1)*x), exp(lambdas(2)*x), exp(lambdas(3)*x)];
P = pinv(X)*y;
P
%P =
%   5.0001 + 0.0000i
%   2.0000 + 0.0001i
%   2.0000 - 0.0001i

Note the cosine multiplier $4$ appears in $P$ as $2$ times $2$.

If we were to fit a $sine$ instead we would get:

y  =  5*exp(0.5*x) + 4*sin(0.15*x);

% other stuff

lambdas =

   0.5000 +      0i
  -0.0000 + 0.1500i
  -0.0000 - 0.1500i

P =

   5.0001e+00 + 5.2915e-16i
  -4.6543e-05 - 1.9999e+00i
  -4.6543e-05 + 1.9999e+00i

Note now the $2$ times $2$ appears in the complex part of $P$.

Exponential multiplied by $x$

If it is desired to fit a sum of exponentials, some with repeated $\lambda$ but multiplied by $x$, again the exact same method can be used. This scenario manifest itself as repeated eigen-values of $\hat{A}$:

clear all;
clc;
% get data
dx = 0.02;
x  = (dx:dx:1.5)';
y  = 5*exp(0.5*x) + 4*exp(-3*x) + 2*x.*exp(-3*x);

% calculate integrals
iy1 = cumtrapz(x, y);
iy2 = cumtrapz(x, iy1);
iy3 = cumtrapz(x, iy2);

% get exponentials lambdas
Y = [iy1, iy2, iy3, x.^2, x, ones(size(x))];
A = pinv(Y)*y;

lambdas = eig([A(1), A(2), A(3); 1, 0, 0; 0, 1, 0]);
lambdas
%lambdas =
%  -2.9991
%  -2.9991 NOTE : repeated means x * e^{lambda * x}
%   0.5000

% get exponentials multipliers
X = [exp(lambdas(1)*x), x.*exp(lambdas(2)*x), exp(lambdas(3)*x)];
P = pinv(X)*y;
P
%P =
%   4.0001
%   1.9955
%   5.0000

What if I don't know how many exponentials fit my data?

If the structure of the model to fit is unknown, one can simply run the algorithm for a large $n$ and compute the singular values of $Y$. By looking at the singular values, we can infer how many of them are actually needed to fit most of the data. From this we can select a minimum $n$ and based in the singular values of $\hat{A}$ we can also deduce the structure of the model (how many exponentials, sines, cosines, if there are exponentials multiplied times $x$, etc.).

clear all;
clc;
% get data
dx = 0.02;
x  = (dx:dx:1.5)';
y  = 5*exp(0.5*x) + 4*exp(-3*x) + 2*exp(-2*x);

% calculate n integrals of y and n-1 powers of x
n = 10;
iy = zeros(length(x), n);
xp = zeros(length(x), n-1);
iy(:,1) = cumtrapz(x, y);
xp(:,1) = ones(size(x));
for ii=2:1:n
    iy(:, ii) = cumtrapz(x, iy(:, ii-1));
    xp(:, ii) = xp(:, ii-1) .* x;
end

% calculate singular values of Y
Y = [iy, xp];
ysv = svd(Y);

% scale singular values to percentage
ysv_scaled = 100 * ysv ./ sum(ysv)
bar(ysv_scaled);
%ysv_scaled =
%   6.9606e+01
%   2.3696e+01
%   4.2812e+00
%   1.5760e+00
%   6.8682e-01
%   1.3106e-01
%   2.0743e-02
%   2.5832e-03
%   2.5537e-04
%   ... even smaller ones

% select n principal components above a threshold of 0.1% 
thres = 0.1;
n = sum(ysv_scaled > 0.1) / 2
% n = 3 (6 singular values) covers about 99.97% of all components contributions
covers = sum(ysv_scaled(1:2*n))
% covers = 99.976

Answer 2

A direct method of fitting (no initial guess, no iteration) for the function : $$y(x)=a+be^{px}+ce^{qx}$$ is summarized below (parameters to be computed : $a,b,c,p,q$ ). It works as well in case of negative $p, q$ :

Instead of minimizing the absolute deviations, the variant below minimizes the relatives deviations :

The theory of this method is given in the paper :https://fr.scribd.com/doc/14674814/Regressions-et-equations-integrales (in French).

The method for the function $$y(x)=a+be^{px}+ce^{qx}+de^{rx}$$ is also available, but not published yet. Contact the author if interrested.

Answer 3

In response to questions raised in comments about the three exponents case, the method of regression based on integral equation is roughly explain below, with a numerical example.

Answer 4

The solution is on wikipedia.

If we want to fit to n exponentials, follow the steps:

Calculate $n$ derivatives of $y$ ; $y^{(n)}, y^{(n-1)}, \cdots , \ddot y , \dot y , y$

Solve linear least squares problem for $[-a_{n-1}, \cdots, -a_2, -a_1, -a_0 ]$ in

$$ y^{(n)} = [y^{(n-1)} , \cdots, \ddot y, \dot y, y] \left[ \begin{array}{c} -a_{n-1} \\ \vdots \\ -a_2 \\ -a_1 \\ -a_0 \end{array} \right] $$

$$ y^{(n)} = Y A \\ Y^T Y A = Y^T y^{(n)} \\ A = (Y^T Y)^{-1}Y^T y^{(n)} $$

Obtain the roots of the characteristic polynomial (using Matlab roots or similar):

$$ z^n + a_{n-1}z^{n-1} + \cdots + a_2 z^2 + a_1z + a_0 $$

The roots $\lambda _n$ are then the values of the exponentials in

$$ y = p_1 e^{\lambda _1 x} + p_2 e^{\lambda _2 x} + \cdots + p_{n-1} e^{\lambda _{n-1} x} + p_n e^{\lambda _n x} $$

To obtain the missing $p_n$, calculate numerically the exponential terms $e^{\lambda _n}$ and solve the least squares problem for $[p_1, p_2, \cdots, p_{n-1}, p_n ]$ in

$$ y = [e^{\lambda _1 x}, e^{\lambda _2 x} , \cdots, e^{\lambda _{n-1} x}, e^{\lambda _n x}] \left[ \begin{array}{c} p_1 \\ p_2 \\ \vdots \\ p_{n-1} \\ p_n \end{array} \right] $$

$$ y = X P \\ X^T X P = X^T y \\ P = (X^T X)^{-1}X^T y $$

Now you have your $n$ exponential model parameters $p_n$ and $\lambda _n$.

I believe this approach is similar of not the same as the one proposed by @JJacquelin, but you have a generic algorithm for $n$ exponentials.

If your $y$ data is noise, I would recommend applying a non-causal low pass filter first.

Here is the Matlab code (you can test it online here):

clear all;
clc;
% get data
dx = 0.01;
x  = (dx:dx:1.5)';
y  =  5*exp(0.5*x) + 4*exp(-3*x) + 2*exp(-2*x);
% calculate derivatives
dy1 = diff(y)/dx;
dy2 = diff(dy1)/dx;
dy3 = diff(dy2)/dx;
% fix derivatives lengths
dy1 = [dy1(1)-(dy1(2)-dy1(1)); dy1];
dy2 = [dy2(1)-2*(dy2(2)-dy2(1)); dy2(1)- (dy2(2)-dy2(1)); dy2];
dy3 = [dy3(1)-3*(dy3(2)-dy3(1)); dy3(1)-2*(dy3(2)-dy3(1)); dy3(1)- (dy3(2)-dy3(1)); dy3];
% get exponentials lambdas
Y = [dy2, dy1, y];
A = pinv(Y)*dy3;
lambdas = roots([1; -A]);
lambdas
%lambdas =
%  -3.0336
%  -1.9933
%   0.4989
% get exponentials multipliers
X = [exp(lambdas(1)*x), exp(lambdas(2)*x), exp(lambdas(3)*x)];
P = pinv(X)*y;
P
%P =
%   3.9461
%   2.0514
%   5.0067

Answer 5

Here is are implementations of the two and three exponent cases, from Jean Jacquelin's work.

#!/usr/bin/python3

'''
Implementation of two and three exponent cases from Jean Jacquelin's
REGRESSION et EQUATION INTEGRALE, available at
https://www.scribd.com/doc/14674814/Regressions-et-equations-integrales
'''    

__author__    = "M Nelson"
__date__      = "January 19, 2020"
__version__   = "0.1"


import sys
import argparse

import numpy as np
from math import sqrt
from scipy.linalg import lstsq
from scipy.optimize import curve_fit

def biexpfit( xdata, ydata, debug=False ):

    x = np.array(xdata)
    y = np.array(ydata)
    if debug:
        print( 'x', x )
        print( 'y', y )

    S = np.empty_like(y)
    S[0] = 0
    S[1:] = np.cumsum(0.5 * (y[1:] + y[:-1]) * np.diff(x))
    if debug:
        print('S', S )

    SS = np.empty_like(y)
    SS[0] = 0
    SS[1:] = np.cumsum(0.5 * (S[1:] + S[:-1]) * np.diff(x))
    if debug:
        print('SS', SS )

    x2 = x * x
    x3 = x2 * x
    x4 = x2 * x2

    M = [ [sum(SS*SS),  sum(SS*S), sum( SS*x2 ), sum(SS*x), sum(SS)],
          [sum(SS*S),   sum(S*S),  sum(S*x2),    sum(S*x), sum(S) ],
          [sum(SS*x2),  sum(S*x2), sum(x4),      sum(x3),  sum(x2) ],
          [sum(SS*x),   sum(S*x),  sum(x3),      sum(x2),  sum(x) ],
          [sum(SS),     sum(S),    sum(x2),      sum(x),   len(xdata) ] ]

    if debug:
        print( 'M' )
        for n in range(5):
            print( M[n] )

    Ycol = np.array( [ sum(SS*y), sum(S*y), sum(x2*y), sum(x*y), sum(y) ] )

    (A,B,C,D,E),residues,rank,singulars = list( lstsq( M, Ycol ) )
    if debug:
        print( 'A-E', A, B, C, D, E )

    '''
    Minv = np.linalg.inv(M)    
    A,B,C,D,E = list( np.matmul(Minv,Ycol) )
    '''

    p = (1/2.)*(B + sqrt(B*B+4*A))
    q = (1/2.)*(B - sqrt(B*B+4*A))
    if debug:
        print( 'p,q', p, q )

    beta = np.exp(p*x)
    eta = np.exp(q*x)

    betaeta = beta * eta

    L = [ [ len(xdata), sum(beta), sum(eta) ],
          [ sum(beta),  sum(beta*beta), sum(betaeta) ],
          [ sum(eta),   sum(betaeta), sum(eta*eta)] ]

    Ycol = np.array( [ sum(y), sum(beta*y), sum(eta*y) ] )

    (a,b,c),residues,rank,singulars = list( lstsq( L, Ycol ) )    
    if debug:
        print( 'a,b,c', a,b,c )

    '''
    Linv = np.linalg.inv(L)
    a,b,c = list( np.matmul( Linv, Ycol ) )
    '''

    # sort in ascending order (fastest negative rate first)
    (b,p),(c,q) = sorted( [[b,p],[c,q]], key=lambda x: x[1])

    return a,b,c,p,q


def triexpfit( xdata, ydata, debug=False ):

    x = np.array(xdata)
    y = np.array(ydata)
    if debug:
        print( 'x', x )
        print( 'y', y )

    S = np.empty_like(y)
    S[0] = 0
    S[1:] = np.cumsum(0.5 * (y[1:] + y[:-1]) * np.diff(x))
    if debug:
        print('S', S )

    SS = np.empty_like(y)
    SS[0] = 0
    SS[1:] = np.cumsum(0.5 * (S[1:] + S[:-1]) * np.diff(x))
    if debug:
        print('SS', SS )

    SSS = np.empty_like(y)
    SSS[0] = 0
    SSS[1:] = np.cumsum(0.5 * (SS[1:] + SS[:-1]) * np.diff(x))
    if debug:
        print('SSS', SSS )

    x2 = x * x
    x3 = x2 * x
    x4 = x3 * x
    x5 = x4 * x
    x6 = x5 * x

    M = [ [ sum(SSS*SSS), sum(SSS*SS), sum(SSS*S), sum(SSS*x3), sum(SSS*x2), sum(SSS*x), sum(SSS) ],
          [ sum(SSS*SS), sum(SS*SS), sum(SS*S), sum(SS*x3), sum(SS*x2), sum(SS*x), sum(SS) ],
          [ sum(SSS*S), sum(SS*S), sum(S*S), sum(S*x3), sum(S*x2), sum(S*x), sum(S) ],
          [ sum(SSS*x3), sum(SS*x3), sum(S*x3), sum(x6), sum(x5), sum(x4), sum(x3) ],
          [ sum(SSS*x2), sum(SS*x2), sum(S*x2), sum(x5), sum(x4), sum(x3), sum(x2) ],
          [ sum(SSS*x), sum(SS*x), sum(S*x), sum(x4), sum(x3), sum(x2), sum(x) ],
          [ sum(SSS), sum(SS), sum(S), sum(x3), sum(x2), sum(x), len(xdata) ] ]

    if debug:
        print( 'M' )
        for n in range(7):
            print( M[n] )

    Ycol = [ sum(y*SSS), sum(y*SS), sum(y*S), sum(y*x3), sum(y*x2), sum(y*x), sum(y) ]

    if debug:
        print( 'Y', Ycol )

    '''
    Minv = np.linalg.inv(M)
    A,B,C,D,E,F,G = list( np.matmul( Minv, Ycol ) )
    print( 'A-G linalg', A, B, C, D, E, F, G )
    '''

    (A,B,C,D,E,F,G),residues,rank,singulars = list( lstsq( M, Ycol ) )
    if debug:
        print( 'A-G', A, B, C, D, E, F, G )

    p,q,r = np.roots( [1.,-C,-B,-A] )
    p,q,r = sorted( [p,q,r] )
    if debug:
        print( 'p,q,r', p, q, r )

    L = [ [ len(xdata), sum(np.exp(p*x)),sum(np.exp(q*x)),sum(np.exp(r*x)) ],
          [ sum(np.exp(p*x)), sum(np.exp(2*p*x)),sum(np.exp((p+q)*x)),sum(np.exp((p+r)*x)) ],
          [ sum(np.exp(q*x)), sum(np.exp((p+q)*x)),sum(np.exp(2*q*x)),sum(np.exp((q+r)*x)) ],
          [ sum(np.exp(r*x)), sum(np.exp((p+r)*x)),sum(np.exp((q+r)*x)),sum(np.exp(2*r*x)) ] ]

    Ycol = [ sum(y), sum(y*np.exp(p*x)),sum(y*np.exp(q*x)),sum(y*np.exp(r*x)) ]

    (a,b,c,d),residues,rank,singulars = list( lstsq( L, Ycol ) )
    if debug:
        print( 'a,b,c,d', a,b,c,d )


    # sort in ascending order (fastest negative rate first)
    (b,p),(c,q),(d,r) = sorted( [[b,p],[c,q],[d,r]], key=lambda x: x[1])

    return a,b,c,d, p,q,r


def __regression_tests__(debug=False):

    x = np.linspace( 0, 2., 100 )

    a, b = 1., 2.
    p = -5.
    print( 'inputs', a, b, p )

    y = a + b * np.exp( p*x )

    a, b, p = expfit( x, y )
    print( 'fits', a, b, p )

    a, b, c = 1.0, .5, 0.25
    p,q = -10., -5.
    print( 'inputs', a,b,c,p,q )

    y = a + b * np.exp( p*x ) + c * np.exp( q* x )

    a,b,c,p,q =  biexpfit( x, y )
    print( 'fits', a,b,c,p,q )


    a, b, c, d = -1., 5., 4., 2.
    p, q, r = .5, -3., -2.
    print( 'inputs', a,b,c,d, p,q,r )

    x = np.linspace( .1, 1.5, 15 )
    print( x )    

    y = a + b * np.exp( p*x ) + c * np.exp( q* x ) + d * np.exp( r* x )
    print( y )

    a,b,c,d, p,q,r = triexpfit( x, y, debug )
    print( 'fits', a,b,c,d, p,q,r )

Answer 6

In order to answer to the comment from Yanqi Huang the method of fitting $$y=be^{px}+ce^{qx}$$ is shown below. This is the same as in the first answer but with one line and column suppressed related to the suppression of the parameter $a$. There is no longer $a\:x^2$ in the linear integral equation.

Answer 7

Fall on your discussion while looking at how to fit Weibull distribution (PDF or CDF). As you may know, these days we look at many exponentials curves which commonly result in the sum of local exponentials. I am no expert but I was giving it a try. I try solving a sum of Weibull distributions (which gives exponentials when the shape parameter equals one). I did not find immediately a minimizing solution, so I was just looking to your method in case. And because a few equations are worth many words:

If any idea or suggestion, gladly appreciated :)

In all cases, thanks so much for this post. I guess I will spend the day implementing your approach.