Fit sum of exponentials 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?
 A: 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.
A: 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.


A: 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

A: 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 )

A: 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.

A: 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. 
