MATLAB function that fits a linear combination of exponential functions to given data

Hello,

I have to write a MATLAB function which looks for a linear combination of exponential functions that fits certain given data the best.

The input of the MATLAB function consists of a (n * 2)-matrix and a vector of initial guesses for the exponents. The data has to fit with the following kind of function: f(x) = C_1*e^(lambda_1 * x) + ... + C_m*e^(lambda_m * x)

So the (n * 2)-matrix has n rows and 2 columns. I suppose the left column contains the x_i and the right column contains the y_i for 1 <= i <= n. The magnitude of m depends on the length of the vector with the initial guesses (lambda_1, ..., lambda_m).

The purpose of this MATLAB function is that the function looks for a function f that is the best choice according to the least squares method. In short, a function such that S = \sum_{i=1}^n (y_i - f(x_i))^2 is minimalized.

The output of the MATLAB function exists of twee vectors with the calculated (best) choices for (lambda_1, ..., lamda_m) and (C_1, ..., C_m) en the residue (thus r_i = y_i - f(x_i)). Furthermore, the function has to make a plot of the data and the fit f.

---

I don't have a problem to realize the output, but I do have trouble with determining f. I can use the function fminsearch, so I can use that for the least squares method I guess.

I did get a hint: For fixed lambda = (lambda_1, ..., lambda_m) I should use the least squares method to find the best coefficients C_j = C_j(lambda).

So I think I have to use the initial guesses of lambda as fixed (at first) and 'only' determine the C's. But I don't have a clue how to do this.

Anyone who can help me? Thanks in advance.

Answer 1

[Image Analyst]

For a fixed lambda, it's really trivial. Just use the standard least squares formula. Did the course cover that? inverse([A'A]) and so on? Just plug in some range of x and get your e^(lambda*x) array and construct your A matrix.

[Tom]

Thanks for your response. I forgot to mention that I am not allowed to use the function lsqnonlin and I think you are talking about that function? If not, what do you mean exactly?

[Image Analyst]

No, I was not thinking of that. This is not a non-linear situation. Look at the equation and the unknowns. If the lambdas are no longer unknown and you have some known, specified values for them, then it turns into a linear situation. The only unknowns in that case are the C's. It doesn't matter that x is plugged into some non-linear equation. The x are not unknowns - they are known. Maybe that's what's tripping you up. The x are your observations - you know those. What matters is if the model equation is linear in the unknowns (which are the C's), which it is. Thus you can use polyfit(). You basically have y = C1*z1 + C2*z2 + .... + Cm*zm. That is an equation which is linear in the C's. You can make up the z vectors from plugging the x's into the exponential formula. Then the z's are columns in the A equation. Go back to the formula for solving a least squares set of equations http://en.wikipedia.org/wiki/Least_squares. The coefficients are solved for by this equation:

C = inverse([z' * z]) * z' * y

Your z is going to be a matrix of the x measurements

z = [e^(lambda_1 * x1), e^(lambda_2 * x1), ... e^(lambda_m * x1);... 
e^(lambda_1 * x2), e^(lambda_2* x2), ... e^(lambda_m * x2);...  
...and so on down until your last x
e^(lambda_1 * xn), e^(lambda_2* xn), ... e^(lambda_m * xn)]

Basically each column of z is the coefficients of the C's if you were to write out all n of your equations in a list. Your instructor should have gone over the derivation of that equation so hopefully this is not the first time you have seen this. (There is a more MATLAB-ish way of solving that like C=y/z if you want to go that route instead of the standard book formula like I showed.)

[Tom]

Thank you very much for your explanation. I now have:

m = length(initialLambda);
datasize = size(data);
n = datasize(1); %number of observations
x = data(:,1);
y = data(:,2);

for i = 1:n,
    for j = 1:m,
        initialz(i,j) = exp(initialLambda(j)*x(i));
    end
end

C = pinv([initialz' * initialz]) * initialz' * y;

f = @(p) 0;
for j = 1:m,
    f = @(p) (f(p) + C(j)*exp(initialLambda(j)*p));
end

for i = 1:n,
    residual(i) = y(i) - f(x(i));
end

optimalLambda = initialLambda;
optimalC = C;

end

This is for the case initialLambda is fixed, so initialLambda == optimalLambda. As output I have optimalLambda, optimalC and residuals. The input consists of data (a (nx2)-matrix of data (x_i y_i)) and the initialLambda (I chose [0.1 0.2 0.3] or something).

=====

Now it works for fixed lambda (the residuals are about 0.22 on average), but I have to modify the MATLAB-function so there isn't a fixed lambda anymore. I guess I have to start with the initialLambda and calculate the C's like I did above and then variate lambda and find the C's again?

I do get as hint: "You can continuously minimize the residue (the residuals in my words) by using fminsearch". I don't understand this. If I look for the smallest residual in this case, I get the x_i and y_i for which the function gives the closest approximation, but that doesn't say a thing about the 'goodness' of the function for all x and y, right?

Thanks in advance.

[Image Analyst]

I haven't used the fminsearch() function. If that's the way they want you do it, then I am not the best person to help you.

[Tom]

It is not necessary to use that approach, but they only gave it as a hint.