Clear Filters
Clear Filters

fminsearch doesn't converge to the right value

20 views (last 30 days)
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:
hold on
[h, ff, s] = sc(0.5, x, y)
plot(x, y)
% 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);
Thanks for help,

Accepted Answer

John D'Errico
John D'Errico on 3 May 2024
Edited: John D'Errico on 3 May 2024
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)';
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!!!!!
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.
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);
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 on 3 May 2024
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
John D'Errico on 3 May 2024
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);
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.

Sign in to comment.

More Answers (1)

Milan Bansal
Milan Bansal on 3 May 2024
Edited: Milan Bansal on 3 May 2024
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);
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.
Hope this helps!
Edoardo_a on 3 May 2024
Thank a lot for the reply! Can't rank your answer caues John replied first...but it was very instructive and helpful!
John D'Errico
John D'Errico on 3 May 2024
Edited: John D'Errico on 3 May 2024
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.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!