algorithm  Matlab  Least Squares data fitting  Cost function with extra constraint
2020腾讯云双十一活动，全年最低！！！（领取3500元代金券），
地址：https://cloud.tencent.com/act/cps/redirect?redirect=1073
【阿里云】双十一活动，全年抄底价，限时3天！(老用户也有)，
入口地址：https://www.aliyun.com/1111/home
I am currently working on some MatLab code to fit experimental data to a sum of exponentials following a method described in this paper.
According to the paper, the data has to follow the following equation (written in pseudocode):
y = sum(v(i)*exp(x/tau(i)),i=1..n)
Here tau(i)
is a set of n
predefined constants. The number of constants also determines the size of the summation, and hence the size of v
. For example, we can try to fit a sum of 100 exponentials, each with a different tau(i)
to our data. However, due to the nature of the fitting and the exponential sum, we need to add another constraint to the problem, and hence to the cost function of the leastsquares method used.
Normally, the cost function of the leastsquares method is given by:
(y_data  sum(v(i)*exp(x/tau(i)),i=1..n)^2
And this has to be minimized. However, to prevent overfitting that would make the timeconstant spectrum extremely noisy, the paper adds the following constraint to the cost function:
v(i)  v(i+1)^2
Because of this extra constraint, as far as I know, the regular algorithms, like lsqcurvefit
aren't useable any longer, and I have to use fminsearch
to search the minimum of my leastsquares cost function with a constraint. The function that has to be minimized, according to me, is the following:
(y_data  sum(v(i)*exp(x/tau(i)),i=1..n)^2 + sum(v(j)  v(j+1)^2,j=1..n1)
My attempt to code this in MatLab is the following. Initially we define the function in a function script, then we use fminsearch
to actually minimize the function and get values for v
.
function res = funcost( v )
%FUNCOST Definition of the function that has to be minimised
%We define a function yvalues with 2 exponentials with known timeconstants
% so we know the result that should be given by minimising.
xvalues = linspace(0,50,10000);
yvalues = 32*exp(xvalues/1)exp(xvalues/10);
%Definition of 30 equidistant point in the logarithmic scale
terms = 30;
termsvector = [1:terms];
tau = termsvector;
for i = 1:terms
tau(i) = 10^(1+3/terms*i);
end
%Definition of the regular function
res_1 = 3;
for i=1:terms
res_1 =res_1+ v(i).*exp(xvalues./tau(i));
end
res_1 = res_1yvalues;
%Added constraint
k=1;
res_2=0;
for i=1:terms1
res_2 = res_2 + (v(i)v(i+1))^2;
end
res=sum(res_1.*res_1) + k*res_2;
end
fminsearch(@funcost,zeros(30,1),optimset('MaxFunEvals',1000000,'MaxIter',1000000))
However, this code is giving me inaccurate results (no error, just inaccurate results), which leads me to believe I either made a mistake in the coding or in the interpretation of the added constraint for the leastsquares method.
algorithm matlab leastsquares datafitting
this question edited Aug 9 '15 at 21:52 Eric Aya 37.6k 14 71 107 asked Aug 9 '15 at 14:31 Spears 6 3 How do you know that the results you get are inaccurate? – A Hernandez Aug 9 '15 at 15:53 The function I use right now to generate the x and yvalues to fit to, 32*exp(t)exp(t/10), should yield increased vamplitudes for tau = 1 and tau = 10. However, this isn't the case when I minimise this function and plot (tau,v). – Spears Aug 9 '15 at 16:01 And that would be the case if you hadn't added the extra term to the cost function. Also fminsearch is a local optimizer, and the solution probably has a local minimum different from the global minimum. Comment out the extra term and give input close to v(i=10)=2 and v(i=20)=1 – A Hernandez Aug 9 '15 at 17:52 Thanks for the suggestion! That might be it (the local minimum different from the global minimum)! If I comment out the extra term, I do get a solution that fits the data. However, then the time domain amplitudes are extremely noisy (around 20004000 for some cases) and switching signs quite often. This is also mentioned in the paper, that the constrain has to be included to prevent this. Because the exponentials don't form an orthogonal basis, there won't be a unique solution I believe, so I think the constrain is needed. – Spears Aug 9 '15 at 18:03

1 Answers
1
I would try to introduce the additional constrain in following way:
res_2 = max((v(1:(end1))v(2:end)).^2);
e.g. instead of minimizing an integrated (summed up) error, it does minmax.
You may also make this constrain stiff by
if res_2 > some_number
k = a_very_big_number;
else
k=0; % or k = a_small_number
end;

this answer edited Aug 10 '15 at 8:02 answered Aug 9 '15 at 17:42 freude 1,330 1 9 28 Thanks for the suggestion! I'll try it out as soon as I can (program is currently running for some different values, so can't try it right away). I'll let you know if it helps! – Spears Aug 9 '15 at 18:04 Tried your constrain! It helps (the code is faster, and I'm getting a visible peak at tau = 10, as expected). However, the peak at 1 is barely visible, which might also be due to a smaller number of points in the range of 10^2..10^0, where this point would be visible. Although the fit is not all that good for lower timeconstants as well (e.g. not going to 0 for t=0), which might be due to the nature of the fitting. – Spears Aug 9 '15 at 22:04
