Optimization Toolbox 4.3

Medium-Scale Nonlinear Data Fitting

This example demonstrates fitting a nonlinear function to data using several of the different medium-scale algorithms available in the Optimization Toolbox™.

Contents

Problem Setup

Consider the following data:

Data = ...
  [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];

Let's plot these data points.

t = Data(:,1);
y = Data(:,2);
axis([0 2 -0.5 6])
hold on
plot(t,y,'ro')
title('Data points')
hold off

We would like to fit the function

   y =  c(1)*exp(-lam(1)*t) + c(2)*exp(-lam(2)*t)

to the data. This function has two linear parameters c and two nonlinear parameters lam.

Since the function has a combination of linear and nonlinear parameters, we will separate the solving into two steps. We will use one of the optimization routines such as LSQNONLIN to solve for the nonlinear parameters, and inside our function we will use "\" to solve for the linear parameters.

We write a function, called FITFUN2, that, given the nonlinear parameters lam and the data, solves for the current estimate of the linear parameters and then returns the error in the fit.

This is the M-file for function FITFUN2:

type fitfun2

function [f, yEst] = fitfun2(lam,Data)
%FITFUN2  Used by DATDEMO to return errors in fitting data to a function.
%   FITFUN2 is used by DATDEMO.
%   f = FITFUN2(lam,Data) returns the error between the data and the values
%   computed by the current function of lam.
%   [f, yEst] = FITFUN2(lam,Data) also returns the estimated value of y;
%   that is, the value of the current model.
%
%   FITFUN2 assumes a function of the form
%
%     y =  c(1)*exp(-lam(1)*t) + ... + c(n)*exp(-lam(n)*t)
%
%   with n linear parameters and n nonlinear parameters.
%
%   To solve for the linear parameters c, we build a matrix A
%   where the j-th column of A is exp(-lam(j)*t) (t is a vector).
%   Then we solve A*c = y for the linear least-squares solution c.

%   Copyright 1990-2006 The MathWorks, Inc.
%   $Revision: 1.1.6.2 $  $Date: 2007/08/03 21:31:49 $

t = Data(:,1); y = Data(:,2);      % separate Data matrix into t and y
A = zeros(length(t),length(lam));  % build A matrix
for j = 1:length(lam)
   A(:,j) = exp(-lam(j)*t);
end
c = A\y;                           % solve A*c = y for linear parameters c
yEst = A*c;
f = y - yEst;                      % compute error (residual) y-A*c

Our objective function requires additional parameters (namely, the matrix Data); the most convenient way to pass these is through an anonymous function:

f = @(x) (norm(fitfun2(x,Data)))

f =

    @(x)(norm(fitfun2(x,Data)))

Next, we set some option parameters via OPTIMSET and we provide a guess for the initial estimates of the nonlinear parameters

options = optimset('LargeScale','off','Display','iter','TolX',1e-3);
lam0 = [1; 0];  % Initial guess for nonlinear parameters

Fit Using Unconstrained Optimization

First, we optimize running the BFGS quasi-Newton algorithm, implemented in the function FMINUNC:

plothandle = plotdatapoints(t,y); % plot data points and get plot handle.
% Output function requires additional parameters data and plothandle; use
% an anonymous function:
foutputfcn = @(x,optimvalues,state) fitfun2outputfcn(x,optimvalues,state, ..
.
                                                     Data,plothandle);
options = optimset(options,'OutputFcn',foutputfcn);
t0 = clock;
[lam,fval,exitflag,output] = fminunc(f,lam0,options);

execution_time = etime(clock, t0);
fprintf('\nNumber of iterations: %g\nNumber of function evaluations: %g\n', 
output.iterations, output.funcCount);
fprintf('Sum of squared residuals at solution: %g\n',fval^2);
fprintf('Execution time: %g\n',execution_time);

                                                        First-order
 Iteration  Func-count       f(x)        Step-size       optimality
     0           3          2.59569                          1.07
     1           6          1.65382       0.936212          0.654
     2           9         0.958437              1          0.266
     3          12           0.8686              1          0.161
     4          18         0.777006       0.227639          0.206
     5          21         0.703335              1          0.368
     6          24         0.537197              1          0.306
     7          27          0.44135              1          0.522
     8          30         0.411956              1          0.487
     9          33         0.386377              1         0.0189
    10          36         0.384929              1         0.0172
    11          39         0.384354              1        0.00133
    12          42         0.384347              1       0.000253
    13          45         0.384347              1      2.02e-005

Local minimum possible.

fminunc stopped because the size of the current step is less than
the selected value of the step size tolerance.




Number of iterations: 13
Number of function evaluations: 45
Sum of squared residuals at solution: 0.147723
Execution time: 0.547

Fit Using Simplex Search

Now we run FMINSEARCH, which implements the Nelder-Mead algorithm:

plothandle = plotdatapoints(t,y); % plot data points and get plot handle.
% Output function requires additional parameters data and plothandle; use
% an anonymous function:
foutputfcn = @(x,optimvalues,state) fitfun2outputfcn(x,optimvalues,state, ..
.
                                                     Data,plothandle);
options = optimset(options,'OutputFcn',foutputfcn);
t0 = clock;
[lam,fval,exitflag,output] = fminsearch(f,lam0,options);
execution_time = etime(clock, t0);
fprintf('\nNumber of iterations: %g\nNumber of function evaluations: %g\n', 
output.iterations, output.funcCount);
fprintf('Sum of squared residuals at solution: %g\n',fval^2);
fprintf('Execution time: %g\n',execution_time);

 Iteration   Func-count     min f(x)         Procedure
     0            1          2.59569
     1            3           2.5425         initial simplex
     2            5          2.51592         expand
     3            7          2.39864         expand
     4            9          2.29017         expand
     5           11          2.02098         expand
     6           13          1.71637         expand
     7           15          1.24546         expand
     8           17          1.03269         expand
     9           18          1.03269         reflect
    10           20          1.03269         contract inside
    11           22          1.03269         contract inside
    12           24          1.03269         contract inside
    13           26          1.03017         contract inside
    14           27          1.03017         reflect
    15           29          1.02989         contract inside
    16           31          1.02989         contract outside
    17           33          1.02981         contract inside
    18           35          1.02974         reflect
    19           37          1.02974         contract inside
    20           39          1.02956         expand
    21           40          1.02956         reflect
    22           42          1.02921         expand
    23           43          1.02921         reflect
    24           45          1.02879         expand
    25           47          1.02817         expand
    26           49          1.02717         expand
    27           51          1.02516         expand
    28           53          1.02283         expand
    29           55          1.01675         expand
    30           57          1.01161         expand
    31           59         0.994259         expand
    32           61         0.979328         expand
    33           63         0.920462         expand
    34           65         0.808253         expand
    35           67         0.550818         expand
    36           69         0.456618         reflect
    37           71         0.456618         contract outside
    38           73         0.456618         contract inside
    39           75         0.456618         contract inside
    40           77         0.456618         contract inside
    41           79         0.454383         reflect
    42           81         0.450431         expand
    43           83         0.438599         expand
    44           85          0.43759         reflect
    45           87         0.423992         expand
    46           89          0.40774         expand
    47           91         0.406488         reflect
    48           93         0.391861         reflect
    49           94         0.391861         reflect
    50           96          0.38652         reflect
    51           98          0.38652         contract inside
    52          100         0.385648         contract outside
    53          102         0.384691         contract inside
    54          104         0.384506         contract inside
    55          106         0.384497         contract outside
    56          108         0.384367         contract inside
    57          110         0.384367         contract inside
    58          112         0.384362         contract outside
    59          114          0.38435         contract inside
    60          116          0.38435         contract inside
    61          118         0.384349         contract outside
    62          120         0.384347         contract inside
    63          122         0.384347         contract inside
    64          123         0.384347         reflect
    65          125         0.384347         contract inside
    66          127         0.384347         contract inside
    67          129         0.384347         contract inside
    68          131         0.384347         contract outside
    69          133         0.384347         contract inside
    70          135         0.384347         contract inside

Optimization terminated:
 the current x satisfies the termination criteria using OPTIONS.TolX of 1.00
0000e-003
 and F(X) satisfies the convergence criteria using OPTIONS.TolFun of 1.00000
0e-004


Number of iterations: 70
Number of function evaluations: 135
Sum of squared residuals at solution: 0.147723
Execution time: 0.593

Fit Using Nonlinear Least Squares (Levenberg-Marquardt)

We now try the Levenberg-Marquardt algorithm in the nonlinear least squares solver LSQNONLIN:

plothandle = plotdatapoints(t,y); % plot data points and get plot handle.
% Output function requires additional parameters data and plothandle; use
% an anonymous function:
F = @(x) fitfun2(x,Data);
foutputfcn = @(x,optimvalues,state) fitfun2outputfcn(x,optimvalues,state, ..
.
                                                     Data,plothandle);
options = optimset(options,'OutputFcn',foutputfcn);
options = optimset(options,'Algorithm','levenberg-marquardt');
t0 = clock;
[lam,resnorm,residual,exitflag,output] = lsqnonlin(F,lam0,[],[],options);
execution_time = etime(clock, t0);
fprintf('\nNumber of iterations: %g\nNumber of function evaluations: %g\n', 
output.iterations, output.funcCount);
fprintf('Sum of squared residuals at solution: %g\n',resnorm);
fprintf('Execution time: %g\n',execution_time);

                                        First-Order                    Norm 
of
 Iteration  Func-count    Residual       optimality      Lambda           st
ep
     0           3         6.73763            2.77         0.01
     1           6         2.10146           0.593        0.001        2.811
18
     2           9         1.28039           0.285       0.0001        1.614
36
     3          15        0.982331           0.292          0.1        1.192
45
     4          18        0.600278           0.147         0.01       0.9572
27
     5          21        0.449997           0.472        0.001        2.362
09
     6          24        0.155291          0.0293       0.0001        2.570
79
     7          27        0.147931          0.0143       1e-005        1.017
04
     8          30        0.147723        0.000195       1e-006       0.1470
52
     9          33        0.147723        3.6e-006       1e-007     0.006346
72

Local minimum possible.

lsqnonlin stopped because the relative size of the current step is less than
the selected value of the step size tolerance.




Number of iterations: 9
Number of function evaluations: 33
Sum of squared residuals at solution: 0.147723
Execution time: 0.297

Fit Using Minimax Optimization

We can also minimize the worst-case error by calling the solver FMINIMAX:

plothandle = plotdatapoints(t,y); % plot data points and get plot handle.
% Output function requires additional parameters data and plothandle; use
% an anonymous function:
F = @(x) fitfun2(x,Data);
foutputfcn = @(x,optimvalues,state) fitfun2outputfcn(x,optimvalues,state, ..
.
                                                     Data,plothandle);
options = optimset(options,'OutputFcn',foutputfcn);
t0 = clock;
options = optimset(options,'MinAbsMax',length(t));

[lam,allfvals,maxfval,exitflag,output]=fminimax(F,lam0,[],[],[],[],[],[],[],
options);
execution_time = etime(clock, t0);
fprintf('\nNumber of iterations: %g\nNumber of function evaluations: %g\n', 
output.iterations, output.funcCount);
fprintf('Sum of squared residuals at solution: %g\n',allfvals'*allfvals);
fprintf('Execution time: %g\n',execution_time);

                  Objective        Max     Line search     Directional
 Iter F-count         value    constraint   steplength      derivative   Pro
cedure
    0      4              0       2.00277
    1      9          1.072        0.1812            1           0.743
    2     14         0.6263        0.1433            1          -0.396
    3     19         0.4369       0.04008            1          -0.189
    4     24         0.3395       -0.0327            1         -0.0755    He
ssian modified
    5     29         0.2947    -0.0001079            1          -0.282    He
ssian modified
    6     34          0.249       0.02101            1          -0.043    He
ssian modified
    7     39         0.2291      0.006745            1         -0.0211
    8     44         0.2272    5.096e-005            1         -0.0301
    9     49         0.2272    5.554e-006            1       -0.000311    He
ssian modified
   10     55         0.2262      0.002457          0.5       -0.000942    He
ssian modified twice
   11     60           0.19       0.01946            1        -0.00698    He
ssian modified
   12     65         0.1791     0.0009595            1         -0.0791    He
ssian modified
   13     70         0.1715     0.0002623            1         -0.0184    He
ssian modified twice
   14     75         0.1715    8.533e-008            1         0.00177    He
ssian modified

Local minimum possible. Constraints satisfied.

fminimax stopped because the size of the current search direction is less th
an
twice the selected value of the step size tolerance and constraints were
satisfied to within the default value of the constraint tolerance.



Active inequalities (to within options.TolCon = 1e-006):
  lower      upper     ineqlin   ineqnonlin
                                     2
                                     5
                                    16

Number of iterations: 15
Number of function evaluations: 75
Sum of squared residuals at solution: 0.199004
Execution time: 2.234