注册 | 登录

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

itPublisher 分享于



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);

%Definition of the regular function
res_1 = 3;

for i=1:terms
    res_1 =res_1+ v(i).*exp(-xvalues./tau(i));

res_1 = res_1-yvalues;

%Added constraint

for i=1:terms-1
    res_2 = res_2 + (v(i)-v(i+1))^2;

res=sum(res_1.*res_1) + k*res_2;



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

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;
    k=0; % or k = a_small_number

  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









您的注册邮箱: 修改

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