Optimization

Function Summary

The following table lists the MATLAB optimization functions.

For more optimization capabilities, see the Optimization Toolbox User's Guide.

Back to Top

Minimizing Functions of One Variable

Given a mathematical function of a single variable coded in an M-file, you can use the fminbnd function to find a local minimizer of the function in a given interval. For example, to find a minimum of the humps.m function in the range (0.3,1), use

x = fminbnd(@humps,0.3,1)

which returns

x =
    0.6370

You can ask for a tabular display of output by passing a fourth argument created by the optimset command to fminbnd

x = fminbnd(@humps,0.3,1,optimset('Display','iter'))

which gives the output

Func-count     x          f(x)         Procedure
    3       0.567376      12.9098        initial
    4       0.732624      13.7746        golden
    5       0.465248      25.1714        golden
    6       0.644416      11.2693        parabolic
    7         0.6413      11.2583        parabolic
    8       0.637618      11.2529        parabolic
    9       0.636985      11.2528        parabolic
   10       0.637019      11.2528        parabolic
   11       0.637052      11.2528        parabolic
 
Optimization terminated:
 the current x satisfies the termination criteria using 
OPTIONS.TolX of 1.000000e-004 


x =

    0.6370

This shows the current value of x and the function value at f(x) each time a function evaluation occurs. For fminbnd, one function evaluation corresponds to one iteration of the algorithm. The last column shows what procedure is being used at each iteration, either a golden section search or a parabolic interpolation.

Back to Top

Minimizing Functions of Several Variables

The fminsearch function is similar to fminbnd except that it handles functions of many variables, and you specify a starting vector x₀ rather than a starting interval. fminsearch attempts to return a vector x that is a local minimizer of the mathematical function near this starting vector.

To try fminsearch, create a function three_var of three variables, x, y, and z.

function b = three_var(v)
x = v(1);
y = v(2);
z = v(3);
b = x.^2 + 2.5*sin(y) - z^2*x^2*y^2;

Now find a minimum for this function using x = -0.6, y = -1.2, and z = 0.135 as the starting values.

v = [-0.6 -1.2 0.135];
a = fminsearch(@three_var,v)

a =
    0.0000   -1.5708    0.1803

Back to Top

fminsearch Algorithm

fminsearch uses the Nelder-Mead simplex algorithm as described in Lagarias et al. [1]. This algorithm uses a simplex of n + 1 points for n-dimensional vectors x. The algorithm first makes a simplex around the initial guess x₀ by adding 5% of each component x₀(i) to x₀, and using these n vectors as elements of the simplex in addition to x₀. (It uses 0.00025 as component i if x₀(i) = 0.) Then, the algorithm modifies the simplex repeatedly according to the following procedure.

Note   The keywords for the fminsearch iterative display appear in bold after the description of the step.

1. Let x(i) denote the list of points in the current simplex, i = 1,...,n+1.

2. Order the points in the simplex from lowest function value f(x(1)) to highest f(x(n+1)). At each step in the iteration, the algorithm discards the current worst point x(n+1), and accepts another point into the simplex. [Or, in the case of step 7 below, it changes all n points with values above f(x(1))].

3. Generate the reflected point

   r = 2m – x(n+1),

   where

   m = Σx(i)/n, i = 1...n,

   and calculate f(r).

4. If f(x(1)) ≤ f(r) < f(x(n)), accept r and terminate this iteration. Reflect

5. If f(r) < f(x(1)), calculate the expansion point s

   s = m + 2(m – x(n+1)),

   and calculate f(s).

   1. If f(s) < f(r), accept s and terminate the iteration. Expand

   2. Otherwise, accept r and terminate the iteration. Reflect

6. If f(r) ≥ f(x(n)), perform a contraction between m and the better of x(n+1) and r:

   1. If f(r) < f(x(n+1)) (i.e., r is better than x(n+1)), calculate

      c = m + (r – m)/2

      and calculate f(c). If f(c) < f(r), accept c and terminate the iteration. Contract outside Otherwise, continue with Step 7 (Shrink).

   2. If f(r) ≥ f(x(n+1)), calculate

      cc = m + (x(n+1) – m)/2

      and calculate f(cc). If f(cc) < f(x(n+1)), accept cc and terminate the iteration. Contract inside Otherwise, continue with Step 7 (Shrink).

7. Calculate the n points

   v(i) = x(1) + (x(i) – x(1))/2

   and calculate f(v(i)), i = 2,...,n+1. The simplex at the next iteration is x(1), v(2),...,v(n+1). Shrink

The following figure shows the points that fminsearch might calculate in the procedure, along with each possible new simplex. The original simplex has a bold outline. The iterations proceed until they meet a stopping criterion.

Back to Top

Example: Curve Fitting

This section gives an example that shows how to fit an exponential function of the form Ae^–λt to some data. The example uses the function fminsearch to minimize the sum of squares of errors between the data and an exponential function Ae^–λt for varying parameters A and λ. This section covers the following topics.

• Creating an M-File for the Example

• Running the Example

• Plotting the Results

Creating an M-File for the Example

To run the example, first create an M-file that:

• Accepts vectors corresponding to the x- and y-coordinates of the data

• Returns the parameters of the exponential function that best fits the data

To do so, copy and paste the following code into an M-file and save it as fitcurvedemo in a directory on the MATLAB path.

function [estimates, model] = fitcurvedemo(xdata, ydata)
% Call fminsearch with a random starting point.
start_point = rand(1, 2);
model = @expfun;
estimates = fminsearch(model, start_point);
% expfun accepts curve parameters as inputs, and outputs sse,
% the sum of squares error for A*exp(-lambda*xdata)-ydata,
% and the FittedCurve. FMINSEARCH only needs sse, but we want
% to plot the FittedCurve at the end.
    function [sse, FittedCurve] = expfun(params)
        A = params(1);
        lambda = params(2);
        FittedCurve = A .* exp(-lambda * xdata);
        ErrorVector = FittedCurve - ydata;
        sse = sum(ErrorVector .^ 2);
    end
end

The M-file calls the function fminsearch, which finds parameters A and lambda that minimize the sum of squares of the differences between the data and the exponential function A*exp(-lambda*t). The nested function expfun computes the sum of squares.

Running the Example

To run the example, first create some random data to fit. The following commands create random data that is approximately exponential with parameters A = 40 and lambda = .5.

xdata = (0:.1:10)'; 
ydata = 40 * exp(-.5 * xdata) + randn(size(xdata));

To fit an exponential function to the data, enter

[estimates, model] = fitcurvedemo(xdata,ydata)

This returns estimates for the parameters A and lambda,

estimates =

   40.1334    0.5025

and a function handle, model, to the function that computes the exponential function A*exp(-lambda*t).

Plotting the Results

To plot the fit and the data, enter the following commands.

plot(xdata, ydata, '*')
hold on
[sse, FittedCurve] = model(estimates);
plot(xdata, FittedCurve, 'r')
 
xlabel('xdata')
ylabel('f(estimates,xdata)')
title(['Fitting to function ', func2str(model)]);
legend('data', ['fit using ', func2str(model)])
hold off

The resulting plot displays the data points and the exponential fit.

Back to Top

Setting Options

You can specify control options that set some minimization parameters using an options structure that you create using the optimset function. You then pass options as an input to the optimization function, for example, by calling fminbnd with the syntax

x = fminbnd(fun,x1,x2,options)

or fminsearch with the syntax

x = fminsearch(fun,x0,options)

For example, to display output from the algorithm at each iteration, set the Display option to 'iter':

options = optimset('Display','iter');

The output structure includes the number of function evaluations, the number of iterations, and the algorithm. The structure appears when you provide fminbnd or fminsearch with a fourth output argument, as in

[x,fval,exitflag,output] = fminbnd(@humps,0.3,1); 

or

[x,fval,exitflag,output] = fminsearch(@three_var,v);

Back to Top

Output Functions

An output function is a function that an optimization function calls at each iteration of its algorithm. Typically, you might use an output function to generate graphical output, record the history of the data the algorithm generates, or halt the algorithm based on the data at the current iteration. You can create an output function as an M-file function, a subfunction, or a nested function.

You can use the OutputFcn option with the following MATLAB optimization functions:

• fminbnd

• fminsearch

• fzero

This section covers the following topics:

• Creating and Using an Output Function

• Structure of the Output Function

• Example of a Nested Output Function

• Fields in optimValues

• States of the Algorithm

• Stop Flag

Creating and Using an Output Function

The following is a simple example of an output function that plots the points generated by an optimization function.

function stop = outfun(x, optimValues, state)
stop = false;
hold on;
plot(x(1),x(2),'.');
drawnow

You can use this output function to plot the points generated by fminsearch in solving the optimization problem

To do so,

1. Create an M-file containing the preceding code and save it as outfun.m in a directory on the MATLAB path.

2. Enter the command

   options = optimset('OutputFcn', @outfun);
   

   to set the value of the Outputfcn field of the options structure to a function handle to outfun.

3. Enter the following commands:

   hold on
   objfun=@(x) exp(x(1))*(4*x(1)^2+2*x(2)^2+x(1)*x(2)+2*x(2));
   [x fval] = fminsearch(objfun, [-1 1], options)
   hold off
   

   This returns the solution

   x =
       0.1290   -0.5323
   
   fval =
      -0.5689
   

   and displays the following plot of the points generated by fminsearch:

   

Structure of the Output Function

The function definition line of the output function has the following form:

stop = outfun(x, optimValues, state)

where

• stop is a flag that is true or false depending on whether the optimization routine should quit or continue. See Stop Flag.

• x is the point computed by the algorithm at the current iteration.

• optimValues is a structure containing data from the current iteration. Fields in optimValues describes the structure in detail.

• state is the current state of the algorithm. States of the Algorithm lists the possible values.

The optimization function passes the values of the input arguments to outfun at each iteration.

Example of a Nested Output Function

The example in Creating and Using an Output Function does not require the output function to preserve data from one iteration to the next. When this is the case, you can write the output function as an M-file and call the optimization function directly from the command line. However, if you want your output function to record data from one iteration to the next, you should write a single M-file that does the following:

• Contains the output function as a nested function—see Nested Functions in MATLAB Programming Fundamentals for more information.

• Calls the optimization function.

In the following example, the M-file also contains the objective function as a subfunction, although you could also write the objective function as a separate M-file or as an anonymous function.

Since the nested function has access to variables in the M-file function that contains it, this method enables the output function to preserve variables from one iteration to the next.

The following example uses an output function to record the points generated by fminsearch in solving the optimization problem

The output function returns the sequence of points as a matrix called history.

To run the example, do the following steps:

1. Open a new M-file in the MATLAB Editor.

2. Copy and paste the following code into the M-file.

   function [x fval history] = myproblem(x0)
       history = [];
       options = optimset('OutputFcn', @myoutput);
       [x fval] = fminsearch(@objfun, x0,options);
           
       function stop = myoutput(x,optimvalues,state);
           stop = false;
           if state == 'iter'
             history = [history; x];
           end
       end
       
       function z = objfun(x)
         z = exp(x(1))*(4*x(1)^2+2*x(2)^2+x(1)*x(2)+2*x(2));
       end
   end
   

3. Save the file as myproblem.m in a directory on the MATLAB path.

4. At the MATLAB prompt, enter

   [x fval history] = myproblem([-1 1])
   

The function fminsearch returns x, the optimal point, and fval, the value of the objective function at x.

x =

    0.1290   -0.5323


fval =

   -0.5689

In addition, the output function myoutput returns the matrix history, which contains the points generated by the algorithm at each iteration, to the MATLAB workspace. The first four rows of history are

history(1:4,:) =

   -1.0000    1.0000
   -1.0000    1.0000
   -1.0750    0.9000
   -1.0125    0.8500 

The final row of points is the same as the optimal point, x.

history(end,:)

ans =

    0.1290   -0.5323

objfun(history(end,:))

ans =

   -0.5689

Fields in optimValues

The following table lists the fields of the optimValues structure that are provided by all three optimization functions, fminbnd, fminsearch, and fzero.

The "Command-Line Display Headings" column of the table lists the headings, corresponding to the optimValues fields that are displayed at the command line when you set the Display parameter of options to 'iter'.

States of the Algorithm

The following table lists the possible values for state:

The following code illustrates how the output function might use the value of state to decide which tasks to perform at the current iteration.

switch state
    case 'init'
          % Setup for plots or guis
    case 'iter'
          % Make updates to plot or guis as needed.
    case 'interrupt'
          % Check conditions to see whether optimization 
          % should quit.
        case 'done'
          % Cleanup of plots, guis, or final plot
otherwise
end

Stop Flag

The output argument stop is a flag that is true or false. The flag tells the optimization function whether the optimization should quit or continue. The following examples show typical ways to use the stop flag.

Stopping an Optimization Based on Data in optimValues.   The output function can stop an optimization at any iteration based on the current data in optimValues. For example, the following code sets stop to true if the objective function value is less than 5:

function stop = myoutput(x, optimValues, state)
stop = false;
% Check if objective function is less than 5.
if optimValues.fval < 5
    stop = true;
end 

Stopping an Optimization Based on GUI Input.   If you design a GUI to perform optimizations, you can make the output function stop an optimization when a user clicks a Stop button on the GUI. The following code shows how to do this, assuming that the Stop button callback stores the value true in the optimstop field of a handles structure called hObject stored in appdata.

function stop = myoutput(x, optimValues, state)
stop = false;
% Check if user has requested to stop the optimization.
stop = getappdata(hObject,'optimstop');

Back to Top

Plot Functions

The PlotFcns field of the options structure specifies one or more functions that an optimization function calls at each iteration to plot various measures of progress while the algorithm executes. Pass a function handle or cell array of function handles. The structure of a plot function is the same as the structure of an output function. For more information on this structure, see Output Functions.

You can use the PlotFcns option with the following MATLAB optimization functions:

• fminbnd

• fminsearch

• fzero

The predefined plot functions for these optimization functions are:

• @optimplotx plots the current point

• @optimplotfval plots the function value

• @optimplotfunccount plots the function count (not available for fzero)

To view or modify a predefined plot function, open the M-file in the MATLAB Editor. For example, to view the M-file for plotting the current point, enter:

edit optimplotx.m

Example: Plot Function

View the progress of a minimization using fminsearch with the plot function @optimplotfval:

1. Write an M-file for the objective function. For this example, use:

   function f = onehump(x)
   
   r = x(1)^2 + x(2)^2;
   s = exp(-r);
   f = x(1)*s+r/20;

2. Set the options to use the plot function:

   options = optimset('PlotFcns',@optimplotfval);

3. Call fminsearch starting from [2,1]:

   [x ffinal] = fminsearch(@onehump,[2,1],options)

4. MATLAB returns the following:

   x =
      -0.6691    0.0000
   
   ffinal =
      -0.4052

   

Back to Top

Finding Roots

The fzero function attempts to find a root of one equation with one variable. You can call this function with either a one-element starting point or a two-element vector that designates a starting interval. If you give fzero a starting point x0, fzero first searches for an interval around this point where the function changes sign. If the interval is found, fzero returns a value near where the function changes sign. If no such interval is found, fzero returns NaN. Alternatively, if you know two points where the function value differs in sign, you can specify this starting interval using a two-element vector; fzero is guaranteed to narrow down the interval and return a value near a sign change.

The following sections contain two examples that illustrate how to find a zero of a function using a starting interval and a starting point. The examples use the function humps.m, which is provided with MATLAB. The following figure shows the graph of humps.

Using a Starting Interval

The graph of humps indicates that the function is negative at x = -1 and positive at x = 1. You can confirm this by calculating humps at these two points.

humps(1)

ans =
    16

humps(-1)

ans =
   -5.1378

Consequently, you can use [-1 1] as a starting interval for fzero.

The iterative algorithm for fzero finds smaller and smaller subintervals of [-1 1]. For each subinterval, the sign of humps differs at the two endpoints. As the endpoints of the subintervals get closer and closer, they converge to zero for humps.

To show the progress of fzero at each iteration, set the Display option to iter using the function optimset.

options = optimset('Display','iter');

Then call fzero as follows:

a = fzero(@humps,[-1 1],options)

This returns the following iterative output:

a = fzero(@humps,[-1 1],options)
 
 Func-count    x          f(x)             Procedure
    2              -1      -5.13779        initial
    3       -0.513876      -4.02235        interpolation
    4       -0.513876      -4.02235        bisection
    5       -0.473635      -3.83767        interpolation
    6       -0.115287      0.414441        bisection
    7       -0.115287      0.414441        interpolation
    8       -0.132562    -0.0226907        interpolation
    9       -0.131666    -0.0011492        interpolation
   10       -0.131618  1.88371e-007        interpolation
   11       -0.131618  -2.7935e-011        interpolation
   12       -0.131618  8.88178e-016        interpolation
   13       -0.131618  8.88178e-016        interpolation
 
Zero found in the interval [-1, 1]

a =

   -0.1316

Each value x represents the best endpoint so far. The Procedure column tells you whether each step of the algorithm uses bisection or interpolation.

You can verify that the function value at a is close to zero by entering

humps(a)

ans =

  8.8818e-016

Using a Starting Point

Suppose you do not know two points at which the function values of humps differ in sign. In that case, you can choose a scalar x0 as the starting point for fzero. fzero first searches for an interval around this point on which the function changes sign. If fzero finds such an interval, it proceeds with the algorithm described in the previous section. If no such interval is found, fzero returns NaN.

For example, if you set the starting point to -0.2, the Display option to Iter, and call fzero by

a = fzero(@humps,-0.2,options)

fzero returns the following output:

Search for an interval around -0.2 containing a sign change:
 Func-count    a          f(a)             b          f(b)      Procedure
    1           -0.2      -1.35385          -0.2     -1.35385   initial interval
    3      -0.194343      -1.26077     -0.205657     -1.44411   search
    5         -0.192      -1.22137        -0.208      -1.4807   search
    7      -0.188686      -1.16477     -0.211314     -1.53167   search
    9         -0.184      -1.08293        -0.216     -1.60224   search
   11      -0.177373     -0.963455     -0.222627     -1.69911   search
   13         -0.168     -0.786636        -0.232     -1.83055   search
   15      -0.154745      -0.51962     -0.245255     -2.00602   search
   17         -0.136     -0.104165        -0.264     -2.23521   search
   18       -0.10949      0.572246        -0.264     -2.23521   search
 
Search for a zero in the interval [-0.10949, -0.264]:
 Func-count    x          f(x)             Procedure
   18        -0.10949      0.572246        initial
   19       -0.140984     -0.219277        interpolation
   20       -0.132259    -0.0154224        interpolation
   21       -0.131617  3.40729e-005        interpolation
   22       -0.131618 -6.79505e-008        interpolation
   23       -0.131618 -2.98428e-013        interpolation
   24       -0.131618  8.88178e-016        interpolation
   25       -0.131618  8.88178e-016        interpolation
 
Zero found in the interval [-0.10949, -0.264]

a =

   -0.1316

The endpoints of the current subinterval at each iteration are listed under the headings a and b, while the corresponding values of humps at the endpoints are listed under f(a) and f(b), respectively.

Note   The endpoints a and b are not listed in any specific order: a can be greater than b or less than b.

For the first nine steps, the sign of humps is negative at both endpoints of the current subinterval, which is shown in the output. At the tenth step, the sign of humps is positive at the endpoint, -0.10949, but negative at the endpoint, -0.264. From this point on, the algorithm continues to narrow down the interval [-0.10949 -0.264], as described in the previous section, until it reaches the value -0.1316.

Back to Top

Troubleshooting

Here is a list of typical problems and recommendations for dealing with them.

options = optimset('FunValCheck', 'on')

Back to Top

Tips

Optimization problems may take many iterations to converge. Most optimization problems benefit from good starting guesses. Providing good starting guesses improves the execution efficiency and may help locate the global minimum instead of a local minimum.

Sophisticated problems are best solved by an evolutionary approach, whereby a problem with a smaller number of independent variables is solved first. Solutions from lower order problems can generally be used as starting points for higher order problems by using an appropriate mapping.

The use of simpler cost functions and less stringent termination criteria in the early stages of an optimization problem can also reduce computation time. Such an approach often produces superior results by avoiding local minima.

Back to Top

Reference

[1] Lagarias, J. C., J. A. Reeds, M. H. Wright, and P. E. Wright. "Convergence Properties of the Nelder-Mead Simplex Method in Low Dimensions." SIAM Journal of Optimization. Vol. 9, Number 1, 1998, pp. 112–147.

Back to Top

Function Handles

Integration

Recommended Products

• MATLAB
• Signal Processing Toolbox
• Optimization Toolbox
• Statistics Toolbox
• Control System Toolbox

Get MATLAB trial software

Includes the most popular MATLAB recorded presentations with Q&A sessions led by MATLAB experts.

Get the Interactive Kit