Documentation Center

  • Trials
  • Product Updates

Unconstrained Nonlinear Optimization Algorithms

Unconstrained Optimization Definition


The term unconstrained means that no restriction is placed on the range of x.

fminunc trust-region Algorithm

Trust-Region Methods for Nonlinear Minimization








where g is the gradient of f at the current point x, H is the Hessian matrix (the symmetric matrix of second derivatives), D is a diagonal scaling matrix, Δ is a positive scalar, and ∥ . ∥ is the 2-norm. Good algorithms exist for solving Equation 6-3 (see [48]); such algorithms typically involve the computation of a full eigensystem and a Newton process applied to the secular equation

Such algorithms provide an accurate solution to Equation 6-3. However, they require time proportional to several factorizations of H. Therefore, for large-scale problems a different approach is needed. Several approximation and heuristic strategies, based on Equation 6-3, have been proposed in the literature ([42] and [50]). The approximation approach followed in Optimization Toolbox solvers is to restrict the trust-region subproblem to a two-dimensional subspace S ([39] and [42]). Once the subspace S has been computed, the work to solve Equation 6-3 is trivial even if full eigenvalue/eigenvector information is needed (since in the subspace, the problem is only two-dimensional). The dominant work has now shifted to the determination of the subspace.



or a direction of negative curvature,




  1. Formulate the two-dimensional trust-region subproblem.

  2. Solve Equation 6-3 to determine the trial step s.


  4. Adjust Δ.



Preconditioned Conjugate Gradient Method


In a minimization context, you can assume that the Hessian matrix H is symmetric. However, H is guaranteed to be positive definite only in the neighborhood of a strong minimizer. Algorithm PCG exits when a direction of negative (or zero) curvature is encountered, i.e., dTHd ≤ 0. The PCG output direction, p, is either a direction of negative curvature or an approximate (tol controls how approximate) solution to the Newton system Hp = –g. In either case p is used to help define the two-dimensional subspace used in the trust-region approach discussed in Trust-Region Methods for Nonlinear Minimization.

fminunc quasi-newton Algorithm

Basics of Unconstrained Optimization

Although a wide spectrum of methods exists for unconstrained optimization, methods can be broadly categorized in terms of the derivative information that is, or is not, used. Search methods that use only function evaluations (e.g., the simplex search of Nelder and Mead [30]) are most suitable for problems that are not smooth or have a number of discontinuities. Gradient methods are generally more efficient when the function to be minimized is continuous in its first derivative. Higher order methods, such as Newton's method, are only really suitable when the second-order information is readily and easily calculated, because calculation of second-order information, using numerical differentiation, is computationally expensive.

Gradient methods use information about the slope of the function to dictate a direction of search where the minimum is thought to lie. The simplest of these is the method of steepest descent in which a search is performed in a direction, –∇f(x), where ∇f(x) is the gradient of the objective function. This method is very inefficient when the function to be minimized has long narrow valleys as, for example, is the case for Rosenbrock's function


The minimum of this function is at x = [1,1], where f(x) = 0. A contour map of this function is shown in the figure below, along with the solution path to the minimum for a steepest descent implementation starting at the point [-1.9,2]. The optimization was terminated after 1000 iterations, still a considerable distance from the minimum. The black areas are where the method is continually zigzagging from one side of the valley to another. Note that toward the center of the plot, a number of larger steps are taken when a point lands exactly at the center of the valley.

Figure 6-1. Steepest Descent Method on Rosenbrock's Function (Equation 6-6)

This function, also known as the banana function, is notorious in unconstrained examples because of the way the curvature bends around the origin. Rosenbrock's function is used throughout this section to illustrate the use of a variety of optimization techniques. The contours have been plotted in exponential increments because of the steepness of the slope surrounding the U-shaped valley.

For an animated version of this figure, enter bandem at the MATLAB® command line.

Quasi-Newton Methods

Of the methods that use gradient information, the most favored are the quasi-Newton methods. These methods build up curvature information at each iteration to formulate a quadratic model problem of the form


where the Hessian matrix, H, is a positive definite symmetric matrix, c is a constant vector, and b is a constant. The optimal solution for this problem occurs when the partial derivatives of x go to zero, i.e.,


The optimal solution point, x*, can be written as


Newton-type methods (as opposed to quasi-Newton methods) calculate H directly and proceed in a direction of descent to locate the minimum after a number of iterations. Calculating H numerically involves a large amount of computation. Quasi-Newton methods avoid this by using the observed behavior of f(x) and ∇f(x) to build up curvature information to make an approximation to H using an appropriate updating technique.

A large number of Hessian updating methods have been developed. However, the formula of Broyden [3], Fletcher [12], Goldfarb [20], and Shanno [37] (BFGS) is thought to be the most effective for use in a general purpose method.

The formula given by BFGS is



As a starting point, H0 can be set to any symmetric positive definite matrix, for example, the identity matrix I. To avoid the inversion of the Hessian H, you can derive an updating method that avoids the direct inversion of H by using a formula that makes an approximation of the inverse Hessian H–1 at each update. A well-known procedure is the DFP formula of Davidon [7], Fletcher, and Powell [14]. This uses the same formula as the BFGS method (Equation 6-10) except that qk is substituted for sk.

The gradient information is either supplied through analytically calculated gradients, or derived by partial derivatives using a numerical differentiation method via finite differences. This involves perturbing each of the design variables, x, in turn and calculating the rate of change in the objective function.

At each major iteration, k, a line search is performed in the direction


The quasi-Newton method is illustrated by the solution path on Rosenbrock's function in Figure 6-2, BFGS Method on Rosenbrock's Function. The method is able to follow the shape of the valley and converges to the minimum after 140 function evaluations using only finite difference gradients.

Figure 6-2. BFGS Method on Rosenbrock's Function

For an animated version of this figure, enter bandem at the MATLAB command line.

Line Search

Line search is a search method that is used as part of a larger optimization algorithm. At each step of the main algorithm, the line-search method searches along the line containing the current point, xk, parallel to the search direction, which is a vector determined by the main algorithm. That is, the method finds the next iterate xk+1 of the form


where xk denotes the current iterate, dk is the search direction, and α* is a scalar step length parameter.

The line search method attempts to decrease the objective function along the line xk + α*dk by repeatedly minimizing polynomial interpolation models of the objective function. The line search procedure has two main steps:

  • The bracketing phase determines the range of points on the line to be searched. The bracket corresponds to an interval specifying the range of values of α.

  • The sectioning step divides the bracket into subintervals, on which the minimum of the objective function is approximated by polynomial interpolation.

The resulting step length α satisfies the Wolfe conditions:


where c1 and c2 are constants with 0 < c1 < c2 < 1.

The first condition (Equation 6-13) requires that αk sufficiently decreases the objective function. The second condition (Equation 6-14) ensures that the step length is not too small. Points that satisfy both conditions (Equation 6-13 and Equation 6-14) are called acceptable points.

The line search method is an implementation of the algorithm described in Section 2-6 of [13]. See also [31] for more information about line search.

Hessian Update

Many of the optimization functions determine the direction of search by updating the Hessian matrix at each iteration, using the BFGS method (Equation 6-10). The function fminunc also provides an option to use the DFP method given in Quasi-Newton Methods (set HessUpdate to 'dfp' in options to select the DFP method). The Hessian, H, is always maintained to be positive definite so that the direction of search, d, is always in a descent direction. This means that for some arbitrarily small step α in the direction d, the objective function decreases in magnitude. You achieve positive definiteness of H by ensuring that H is initialized to be positive definite and thereafter (from Equation 6-15) is always positive. The term is a product of the line search step length parameter αk and a combination of the search direction d with past and present gradient evaluations,


You always achieve the condition that is positive by performing a sufficiently accurate line search. This is because the search direction, d, is a descent direction, so that αk and the negative gradient –∇f(xk)Td are always positive. Thus, the possible negative term –∇f(xk+1)Td can be made as small in magnitude as required by increasing the accuracy of the line search.

fminsearch Algorithm

fminsearch uses the Nelder-Mead simplex algorithm as described in Lagarias et al. [57]. This algorithm uses a simplex of n + 1 points for n-dimensional vectors x. The algorithm first makes a simplex around the initial guess x0 by adding 5% of each component x0(i) to x0, and using these n vectors as elements of the simplex in addition to x0. (It uses 0.00025 as component i if x0(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 = 2mx(n+1),


    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(mx(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 + (rm)/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.

Was this topic helpful?