Fit sum of exponentials

Solution 1:

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.

Solution 2:

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

Solution 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.

Solution 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