Fit sum of multivariate exponentials I am working on a curve fitting task for a function which is multivariate sum of two exponentials.
The function which I am trying to fit has below form:
$$z(x,y)=ae^{bx}y+ce^{dx}$$
and want to find a, b, c, d values given values of x, y and z from measurement data.
Has anyone already worked on such function?
It seems not possible to linearize the problem using function integrals as done for fitting of sum of exponentials with single variable because second variable y is present in only one exponential term.
 A: I came up with a solution exactly at 4:20. Your specific case is a lucky one because the term that depends on the independent variable $y$ affects linearly the non-linear exponentials that depend on $x$. So the problem remains linear and you can still apply the method I proposed in this other answer (with some minor changes).
First you have to store your data $z(x,y)$ in a (tall) column vector, putting each $y$ sample head-to-tail with each other. Assuming you have $nx$ number of samples in $x$, and $ny$ number of samples in $y$:
$$
z = \left[
\begin{array}{c}
z(x, y_1) \\
z(x, y_2) \\
\vdots \\
z(x, y_{ny})
\end{array}
\right]
$$
Now we need to equate the right hand side of the equation according to the proposed method. For all $ny$ samples in $y$, the exponental terms ("lambdas") are the same, but the terms that multiply $x$ and the constant are not, so the matrix should have the form:
$$
Z = \left[
\begin{array}{c}
\int{z(x, y_1)} & \int^2{z(x, y_1)} & x & 1 & 0 & 0 & \dots & 0 & 0 \\
\int{z(x, y_2)} & \int^2{z(x, y_2)} & 0 & 0 & x & 1 & \dots & 0 & 0 \\
\vdots \\
\int{z(x, y_{ny})} & \int^2{z(x, y_{ny})} & 0 & 0 & 0 & 0 & \dots & x & 1
\end{array}
\right]
$$
So now we have:
$$
z = Z A
$$
So:
$$
A = (Z^T Z)^{-1}Z^T z 
$$
$$
\lambda = eig\left(
\left[
\begin{array}{c}
A(1) & A(2) \\
1 & 0
\end{array}
\right]
\right) 
$$
The rest follows by the proposed method. The respective Matlab code is:
clear all;
clc;

%% Create test data

% independent variable x
dx = 0.02;
x  = (dx:dx:1.5)';
nx = length(x);
% independent variable y
dy = 0.05;
y  = (dy:dy:1.5)';
ny = length(y);
% dependent variable z(x, y)
% stored in column vector head-to-tail for each y[k]
z = zeros(nx*ny, 1);
for k = 1:1:ny
    z_start = nx*(k-1)+1;
    z_end   = nx*k;
    % calculate z(y[k])
    z(z_start:z_end) = 5*exp(0.5*x)*y(k) + 4*exp(-3*x);
end

%% Estimate exponential contants ("lambdas")

% build least squares matrix Z to obtain lambdas
Z = zeros(nx*ny, 2+2*ny);
for k = 1:1:ny
    z_start = nx*(k-1)+1;
    z_end   = nx*k;    
    % calculate integrals of z(y[k]) wrt x, for each kth sample in y
    iz1 = cumtrapz(x, z(z_start:z_end));
    iz2 = cumtrapz(x, iz1);
    % store in appropriate location
    Z(z_start:z_end, 1:2) = [iz1, iz2];
    Z(z_start:z_end, 2+2*(k-1)+1:2+2*(k-1)+2) = [x, ones(size(x))];
end
% get exponentials lambdas
A = pinv(Z)*z;
lambdas = eig([A(1), A(2); 1, 0]);
lambdas
%lambdas =
%  -2.9991
%   0.5000

%z_hat = Z*A;
%figure();
%plot(z, 'b-'); hold on;
%plot(z_hat, 'r--'); hold on;

%% Estimate exponentials multipliers

% build least squares matrix X to obtain multipliers
X = zeros(nx*ny, 2);
for k = 1:1:ny
    x_start = nx*(k-1)+1;
    x_end   = nx*k;    
    X(x_start:x_end, :) = [exp(lambdas(2)*x)*y(k), exp(lambdas(1)*x)];
end
P = pinv(X)*z;
P
%P =
% 4.9999
% 3.9995

You can try it online here:
https://octave-online.net/bucket~UfZSeFmySresJFMLNy4Y34
Please upvote my proposed method in the other question if you found it useful.

Edit:
@PayasVartak, the lambdas computation using the eigenvalues is indeed the optimal solution for the two lambdas (there only two values). The reason why you you might not be getting the best fit, might be because of the second least squares problem, in the "Estimate exponentials multipliers" section of the code.
Note that I used the equation:
[exp(lambdas(2)*x)*y(k), exp(lambdas(1)*x)]

for the fit, because I already knew in advance which exponential gets multiplied by y(k). But actually we don't know at this point which of the two exp(lambda(?)*x) is the one that gets multiplied by y(k). So actually you need to try also:
[exp(lambdas(2)*x), exp(lambdas(1)*x)*y(k)]

And see which one makes the best fit.
As @JJacquelin proposed, it might be a good idea to post some of your experimental data so we can help you better.
A: The method of using integral equation to linearize the problem is explained in https://fr.scribd.com/doc/14674814/Regressions-et-equations-integrales . The case of two exponentials is treated pp.71-74. More exponentials are possible as well.
But the theory is done only in case of one variable.
In this paper the case of two variables is passed over superficially.pp.75-80. It is shown that the method can be extended in case of two variables if the data has a specific structure as shown in the numerical example pp.81-84 .
The case of $z(x,y)=ae^{bx}y+ce^{dx}$ is not explicitly treated. This would need more studies. The structure of the data (How the experimental points are presented in a table) is very important in practical application. In order to give a preliminary  advice on the feasibility in your case, you should joint to your question an example of data.
