ITKeyword,专注技术干货聚合推荐

注册 | 登录

algorithm - Matlab - Least Squares data fitting - Cost function with extra constraint

itPublisher 分享于

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 pseudo-code):

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 least-squares method used.

Normally, the cost function of the least-squares 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 over-fitting that would make the time-constant 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 least-squares 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..n-1)

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 time-constants
% so we know the result that should be given by minimising.
xvalues = linspace(0,50,10000);
yvalues = 3-2*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_1-yvalues;

%Added constraint
k=1;
res_2=0;

for i=1:terms-1
    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 least-squares method.

algorithm matlab least-squares data-fitting
|
  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 y-values to fit to, 3-2*exp(-t)-exp(-t/10), should yield increased v-amplitudes 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 2000-4000 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:(end-1))-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 time-constants 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

 | 


相关阅读排行


相关内容推荐

最新文章

×

×

请激活账号

为了能正常使用评论、编辑功能及以后陆续为用户提供的其他产品,请激活账号。

您的注册邮箱: 修改

重新发送激活邮件 进入我的邮箱

如果您没有收到激活邮件,请注意检查垃圾箱。