Nonlinear Data-Fitting
This example shows how to fit a nonlinear function to data using several Optimization Toolbox™ algorithms.
Problem Setup
The problem is to fit a function to the following data.
xdata = ...
[0.0000 5.8955
0.1000 3.5639
0.2000 2.5173
0.3000 1.9790
0.4000 1.8990
0.5000 1.3938
0.6000 1.1359
0.7000 1.0096
0.8000 1.0343
0.9000 0.8435
1.0000 0.6856
1.1000 0.6100
1.2000 0.5392
1.3000 0.3946
1.4000 0.3903
1.5000 0.5474
1.6000 0.3459
1.7000 0.1370
1.8000 0.2211
1.9000 0.1704
2.0000 0.2636];
Plot the data points.
t = xdata(:,1); y = xdata(:,2); plot(t,y,"o") title("Data points")
Try to fit the function
to the data.
Solution Approach Using lsqcurvefit
To fit a function to data using lsqcurvefit
, map the parameters in the function to a variable x
as follows:
x(1) =
x(2) =
x(3) =
x(4) =
Define the parameters to fit, x
, in a function handle along with the data xdata
.
F = @(x,xdata)x(1)*exp(-x(2)*xdata) + x(3)*exp(-x(4)*xdata);
Arbitrarily set the initial point x0 as follows: = 1, = 1, = 1, = 0.
x0 = [1 1 1 0];
Run the lsqcurvefit
solver and plot the resulting fit.
[x,resnorm,~,exitflag,output] = lsqcurvefit(F,x0,t,y)
Local minimum possible. lsqcurvefit stopped because the final change in the sum of squares relative to its initial value is less than the value of the function tolerance. <stopping criteria details>
x = 1×4
3.0068 10.5869 2.8891 1.4003
resnorm = 0.1477
exitflag = 3
output = struct with fields:
firstorderopt: 7.8830e-06
iterations: 6
funcCount: 35
cgiterations: 0
algorithm: 'trust-region-reflective'
stepsize: 0.0096
message: 'Local minimum possible.↵↵lsqcurvefit stopped because the final change in the sum of squares relative to ↵its initial value is less than the value of the function tolerance.↵↵<stopping criteria details>↵↵Optimization stopped because the relative sum of squares (r) is changing↵by less than options.FunctionTolerance = 1.000000e-06.'
bestfeasible: []
constrviolation: []
hold on plot(t,F(x,t)) hold off
Solution Approach Using fminunc
To solve the problem using fminunc
, define the objective function as the sum of squares of the residuals.
Fsumsquares = @(x)sum((F(x,t) - y).^2);
[xunc,ressquared,eflag,outputu] = ...
fminunc(Fsumsquares,x0)
Local minimum found. Optimization completed because the size of the gradient is less than the value of the optimality tolerance. <stopping criteria details>
xunc = 1×4
2.8890 1.4003 3.0069 10.5862
ressquared = 0.1477
eflag = 1
outputu = struct with fields:
iterations: 30
funcCount: 185
stepsize: 0.0017
lssteplength: 1
firstorderopt: 2.9662e-05
algorithm: 'quasi-newton'
message: 'Local minimum found.↵↵Optimization completed because the size of the gradient is less than↵the value of the optimality tolerance.↵↵<stopping criteria details>↵↵Optimization completed: The first-order optimality measure, 9.409720e-07, is less ↵than options.OptimalityTolerance = 1.000000e-06.'
Notice that fminunc
finds the same solution as lsqcurvefit
, but takes many more function evaluations to do so. The parameters for fminunc
are in the opposite order as those for lsqcurvefit
; the larger is , not . This should not be surprising, because the order of variables is arbitrary.
fprintf(['There were %d iterations using fminunc,' ... ' and %d using lsqcurvefit.\n'], ... outputu.iterations,output.iterations)
There were 30 iterations using fminunc, and 6 using lsqcurvefit.
fprintf(['There were %d function evaluations using fminunc,' ... ' and %d using lsqcurvefit.'], ... outputu.funcCount,output.funcCount)
There were 185 function evaluations using fminunc, and 35 using lsqcurvefit.
Splitting the Linear and Nonlinear Problems
Notice that the fitting problem is linear in the parameters and . This means for any values of and , you can use the backslash operator to find the values of and that solve the least-squares problem.
Rework the problem as a two-dimensional problem, searching for the best values of and . The values of and are calculated at each step using the backslash operator as described above.
function yEst = fitvector(lam,xdata,ydata) %FITVECTOR Used by DATDEMO to return value of fitting function. % yEst = FITVECTOR(lam,xdata) returns the value of the fitting function, y % (defined below), at the data points xdata with parameters set to lam. % yEst is returned as a N-by-1 column vector, where N is the number of % data points. % % FITVECTOR assumes the fitting function, y, takes the form % % y = c(1)*exp(-lam(1)*t) + ... + c(n)*exp(-lam(n)*t) % % with n linear parameters c, and n nonlinear parameters lam. % % To solve for the linear parameters c, we build a matrix A % where the j-th column of A is exp(-lam(j)*xdata) (xdata is a vector). % Then we solve A*c = ydata for the linear least-squares solution c, % where ydata is the observed values of y. A = zeros(length(xdata),length(lam)); % Build the A matrix. for j = 1:length(lam) A(:,j) = exp(-lam(j)*xdata); end c = A\ydata; % Solve A*c = y for linear parameters c. yEst = A*c; % Return the estimated response based on c. end
Solve the problem using lsqcurvefit
, starting from a two-dimensional initial point []:
x02 = [1 0]; F2 = @(x,t) fitvector(x,t,y); [x2,resnorm2,~,exitflag2,output2] = lsqcurvefit(F2,x02,t,y)
Local minimum possible. lsqcurvefit stopped because the final change in the sum of squares relative to its initial value is less than the value of the function tolerance. <stopping criteria details>
x2 = 1×2
10.5861 1.4003
resnorm2 = 0.1477
exitflag2 = 3
output2 = struct with fields:
firstorderopt: 4.4087e-06
iterations: 10
funcCount: 33
cgiterations: 0
algorithm: 'trust-region-reflective'
stepsize: 0.0080
message: 'Local minimum possible.↵↵lsqcurvefit stopped because the final change in the sum of squares relative to ↵its initial value is less than the value of the function tolerance.↵↵<stopping criteria details>↵↵Optimization stopped because the relative sum of squares (r) is changing↵by less than options.FunctionTolerance = 1.000000e-06.'
bestfeasible: []
constrviolation: []
The efficiency of the two-dimensional solution is similar to that of the four-dimensional solution:
fprintf(['There were %d function evaluations using the 2-d ' ... 'formulation, and %d using the 4-d formulation.'], ... output2.funcCount,output.funcCount)
There were 33 function evaluations using the 2-d formulation, and 35 using the 4-d formulation.
Split Problem is More Robust to Initial Guess
Choosing a bad starting point for the original four-parameter problem leads to a local solution that is not global. However, choosing a starting point with the same bad and values for the split two-parameter problem leads to the global solution. To demonstrate this result, re-run the original problem with a start point that leads to a relatively bad local solution, and compare the resulting fit with the global solution.
x0bad = [5 1 1 0];
[xbad,resnormbad,~,exitflagbad,outputbad] = ...
lsqcurvefit(F,x0bad,t,y)
Local minimum possible. lsqcurvefit stopped because the final change in the sum of squares relative to its initial value is less than the value of the function tolerance. <stopping criteria details>
xbad = 1×4
-22.9036 2.4792 28.0273 2.4791
resnormbad = 2.2173
exitflagbad = 3
outputbad = struct with fields:
firstorderopt: 0.0056
iterations: 32
funcCount: 165
cgiterations: 0
algorithm: 'trust-region-reflective'
stepsize: 0.0021
message: 'Local minimum possible.↵↵lsqcurvefit stopped because the final change in the sum of squares relative to ↵its initial value is less than the value of the function tolerance.↵↵<stopping criteria details>↵↵Optimization stopped because the relative sum of squares (r) is changing↵by less than options.FunctionTolerance = 1.000000e-06.'
bestfeasible: []
constrviolation: []
hold on plot(t,F(xbad,t),'g') legend("Data","Global fit","Bad local fit","Location","NE") hold off
fprintf(['The residual norm at the good ending point is %f,\r' ... 'and the residual norm at the bad ending point is %f.'], ... resnorm,resnormbad)
The residual norm at the good ending point is 0.147723, and the residual norm at the bad ending point is 2.217300.