fminsearch doesn't converge to the right value

Hi all,

I am trying to solve and minimize the following exponential fit, trying to retrieve the exponential term.

I am using fminsearch to minimize the "sum-squared-error cost function" but somehow it is not really working. I am sure I am making some stupid mistake but I can't figure it out ... It should give me back the 0.1 value but never really get it (even if noise is 0).

Here the code:

x=0:100;
y=exp(-0.1*x)'+0.1*rand(101,1);
plot(x,y)
hold on
%%
[h, ff, s] = sc(0.5, x, y)
plot(x, y)
plot(x,h,x,y,x,ff)
% hold on
% s
function [h ff s]=sc(k, x, y)
co=[exp(-k*x)' ones(size(y))]\y;
%
% ff=co(1)*exp(-k*x)'+co(2)*ones(size(y));
% h=y-ff;
% s=sum(h.^2);
%
ff = @(k,x)  co(1).*exp(-k.*x)'+co(2).*ones(size(y));
h  = @(k) sum((y-ff(k,x)).^2);
options = optimset('PlotFcns',@optimplotfval,'Display','iter');
s  = fminsearch(h,k,options);
ff=exp(-s*x)'+co(2)*ones(size(y));
h=y-ff;
end

Thanks for help,

E

Answer 1

[John D'Errico]

These things are always easier than you think. Well, except when they are just as hard as you think. ;-)

First, I'll create x as a COLUMN vector. That will avoid problems later, since you were forced to constantly transpose it. Just start out right from the beginning!

x = (0:100)';

Next, WHY IN THE NAME OF GOD AND LITTLE GREEN APPLES DID YOU USE RAND??????????????????

Do you understand that rand creates random numbers in the interval [0,1]? There will NEVER BE NEGATIVE NUMBERS COMING FROM RAND. USE RANDN INSTEAD!

y = exp(-0.1*x) + 0.1*rand(101,1); % THIS IS JUST BAD!!!!!
plot(x,y)

Next, plot it. Look at what you have. Note that here, the negative exponential is pretty much garbage noise for most of the curve, so the only useful signal comes from the initial part. But, ok. You still have enough data to get a fit out of it.

The model you have chosen is of the general form

y = a*exp(-k*x) + b

You MIGHT expect to find a==1, b == 0, and k = 0.1. Except that you used rand! And that means I will expect to see approximately b = 0.05.

Of course, those parameters will be estimated from somewhat poor data, but we will still get a valid result. Just for kicks, I'l throw the curve fitting toolbox at it first.

mdl = fittype('a*exp(-k*x) + b','indep','x')

mdl = 
     General model:
     mdl(a,b,k,x) = a*exp(-k*x) + b

fittedmdl = fit(x,y,mdl)

Warning: Start point not provided, choosing random start point.

fittedmdl = 
     General model:
     fittedmdl(x) = a*exp(-k*x) + b
     Coefficients (with 95% confidence bounds):
       a =       1.016  (0.9837, 1.048)
       b =     0.06034  (0.05328, 0.0674)
       k =      0.1028  (0.0972, 0.1083)

Even without starting values provided, fit was happy. We got something reasonable. Note that the estimate of b would be expected to be roughly 0.05. NOT zero. This is because the mean of a random number with a uniform distribution on the interval [0.0.1] is 0.05.

plot(fittedmdl,x,y)

I would not expect anything better. Now to look at your use of fminsearch for the fit. It appears you are trying to use a partitioned least squares, where you separate the conditionally linear parameters from the linear parameters. But you don't understand how to write that. (You could just download my fminspleas from the file exchange.)

k0 = 0.25;
[kest,abest]= pleasFit(k0, x, y)

 
 Iteration   Func-count         f(x)         Procedure
     0            1         0.919486         
     1            2         0.919486         initial simplex
     2            4         0.830404         expand
     3            6         0.611001         expand
     4            8         0.343807         reflect
     5           10         0.278758         contract outside
     6           12         0.278758         contract inside
     7           14         0.278758         contract inside
     8           16         0.277347         contract inside
     9           18         0.277347         contract inside
    10           20         0.277347         contract inside
    11           22         0.277324         contract inside
    12           24         0.277324         contract inside
    13           26         0.277324         contract inside
 
Optimization terminated:
 the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-04 
 and F(X) satisfies the convergence criteria using OPTIONS.TolFun of 1.000000e-04 

kest = 0.1027

abest = 2x1

1.0159 0.0603

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

As you can see, it returns the same estimates as did fit. Again, the estimate of b will be screwy because of your use of rand. The estimate I would espect to see for this model, with this data is ....

y = 1*exp(-0.1*x) + 0.05

function [kest,abest] = pleasFit(k0, x, y)
% Partitioned least squares estimation for the model
%   y = a*exp(-k*x) + b
%
% k will be provided as an initial value, x and y are column vectors
linearest = @(x,y,k) [exp(-k*x),ones(size(y))]\y;
obj = @(k) norm([exp(-k*x),ones(size(y))]*linearest(x,y,k) - y);
options = optimset('PlotFcns',@optimplotfval,'Display','iter');
[kest,err,exitflag] = fminsearch(obj,k0,options);
abest = linearest(x,y,kest);
end

In the end, there were two MAJOR issues. First, WHY DID YOU USE RAND? RAND generates UNIFORM numbers in the interval [0,1]. I think too often, people just assume a random number is a random number, and do not recognize the difference between them. The result is you will get confusing results, because the additive random part does not have a mean of zero when you used rand.

Next, your code for the fit using fminsearch was completely wrong. It needs to estimate the parameters a and b INSIDE the objective. So for every value of k that fminsearch tries, it needs to estimate a NEW set of estimates for a and b. This was your major problem of course.

A partitioned nonlinear least squares splits the problem into two sets of unknowns. The nonlinear parameter (k), and then two conditionally linear parameters, here a and b. If you knew the value of k, you could estimate a and b. In fact, you could then just use polyfit to estimate a and b. (I probably should have done so.) But the linear part of the fit needs to happen INSIDE the objective function.

[Edoardo_a]

Hi John, thanks a lot for the deep explanation, I am still learning how to use fminsearch ...and sorry for using "rand", lesson learnt!

You mean uinsg polyfit instead of linearest? Ok...but this work very well too, Thanks a lot!

[John D'Errico]

Exactly. I could have written it like this instead:

function [kest,abest] = pleasFit(k0, x, y)
% Partitioned least squares estimation for the model
%   y = a*exp(-k*x) + b
%
% k will be provided as an initial value, x and y are column vectors
obj = @(k) norm(polyval(polyfit(exp(-k*x),y,1),exp(-k*x)) - y);
options = optimset('PlotFcns',@optimplotfval,'Display','iter');
[kest,err,exitflag] = fminsearch(obj,k0,options);
abest = polyfit(exp(-k*x),y,1) - y);
end

I'll concede that few people even attempt ideas like partitioned least squares. So kudos to you for making a good stab at it.

Once you get the hang of it, the process works very well. Partitioned least squares is far more robust than a general nonlinear least squares, since it greatly reduces the size of the parameter space. Here, we have only a 1 variable problem to search using fminsearch. That makes the search considerably faster, since you are only searching over a 1-d parameter space. It make the problem more robust, since there is no need to provide 3 starting values.

Answer 2

[Milan Bansal]

I understand that you are facing some issues when using "fminsearch" for optimizing the "sum-squared-error" cost function to find the best fit for the exponential model , where A is a constant, B is noise, and "k" is the exponential decay rate you are trying to estimate.

Please refer to the following points to estimate the value of "k":

• Define an objective function that calculates the sum-squared-error between the model predictions and the actual data for a given "k". This function will be minimized to find the best "k".
• Use the fminsearch function to find the "k" that minimizes the objective function. Start with an initial guess for "k" and let fminsearch adjust "k" to minimize the sum-squared-error.

Here is how you can modify your code to resolve the issue:.

x = 0:100;
y = exp(-0.1*x)' + 0.1*randn(101,1); 
plot(x, y, 'b.'); % Plotting the original data
hold on;

initialGuess = 0.5; % Initial guess for k
[kOpt, fval] = optimizeK(initialGuess, x, y); % Optimization to find the best k

 
 Iteration   Func-count         f(x)         Procedure
     0            1          4.31967         
     1            2          4.31967         initial simplex
     2            4          4.08769         expand
     3            6          3.49604         expand
     4            8          1.51058         expand
     5           10          1.51058         contract inside
     6           12          1.04957         contract outside
     7           14          1.04957         contract inside
     8           16          1.04957         contract outside
     9           18          1.04254         contract inside
    10           20          1.04254         contract inside
    11           22          1.04218         contract inside
    12           24          1.04217         contract inside
    13           26          1.04213         contract inside
    14           28          1.04213         contract inside
    15           30          1.04213         contract inside
 
Optimization terminated:
 the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-04 
 and F(X) satisfies the convergence criteria using OPTIONS.TolFun of 1.000000e-04 

% Generating the fitted curve with the optimized k
yFitted = exp(-kOpt*x);
plot(x, yFitted, 'r-'); % Plotting the fitted curve
title(sprintf('Fitted Curve with k = %f', kOpt));
hold off;

function [kOpt, fval] = optimizeK(kInitial, x, y)
    objectiveFunction = @(k) sum((y - exp(-k*x)').^2); % Objective function to minimize sum-squared-error

    % Using fminsearch to find the optimal k
    options = optimset('PlotFcns',@optimplotfval, 'Display', 'iter');
    [kOpt, fval] = fminsearch(objectiveFunction, kInitial, options);
end

The value of "k" estimated from above code is 0.094, which is close to the desired value 0.1. The difference is due to the added noise. If the noise become zero, resultant "k" will be 0.1.

Please refer to the following documentation to learn more about fminsearch.

https://www.mathworks.com/help/matlab/ref/fminsearch.html

Hope this helps!

[Edoardo_a]

Thank a lot for the reply! Can't rank your answer caues John replied first...but it was very instructive and helpful!

Thanks!

[John D'Errico]

Argh. NO!!!! This answer was completely wrong! It did not use a partitioned least squares estimator! It just estimated k, given the model

y= exp(-k*x)

You need to see the difference. What was done here was @Milan Bansal did nothing more than assume the above model, even though initially it ws stated the model would be

y = A*exp(-k*x) + B

I did say before that very few people seem to appreciate what a partitioned least squares estimation is and how it works. This is a good example.