Technical Articles

Optimally Speaking

By Cleve Moler, MathWorks


Version 2.0 of the Optimization Toolbox, included on Release 11 of the MathWorks CD, is a major upgrade, featuring new methods for large-scale, sparse optimization problems. The original version of the toolbox, written several years ago by Andy Grace, has proved useful for the small- and medium-scale problems often encountered in optimal control. Version 2.0 extends the capabilities in a number of ways. Tom Coleman of Cornell University and Mary Ann Branch of The MathWorks have added new functions for nonlinear optimization problems involving hundreds or thousands of variables and constraints. And Penny Anderson, also of The MathWorks, has adapted the LIPSOL routine to solve large-scale linear programming problems. The LIPSOL routine was originally written by Yin Zhang, who is now at Rice University.

For many years, large-scale optimization problems were limited to linear programs, which are models with linear objective functions and linear constraints. The simplex algorithm, developed by George Dantzig in the 1940s, proved to be much more effective in practice than its theoretical worst case computational complexity would suggest. In 1984, Narendra Karmarkar described an algorithm with much better theoretical complexity properties. Karmarkar also made strong claims about the effectiveness of his algorithm but refused to disclose the details to peer review. Other researchers were eventually able to relate Karmarkar's method to more classic methods, improve upon it, and develop publicly accessible codes.

For linear problems involving n variables, the feasible region is the polygonal subset of n-dimensional space where the equality and inequality constraints are satisfied. The optimizing solution must lie on the boundary of the feasible region, so the simplex algorithm limits its search to this boundary. Karmarkar's innovation was an algorithm that searches over the interior of the feasible region and only approaches the boundary as the iteration converges.

LIPSOL is Zhang's MATLAB implementation of the linear programming techniques that have resulted from the research on interior point methods. The algorithm is now described as the Mehrotra predictor-corrector variant of a primal-dual interior-point method. The implementation makes use of MATLAB's built-in sparse matrix operations. For large linear programming problems, this new linprog function may be orders of magnitude more efficient than it was in the first version of the toolbox.

Recent years have also seen significant improvements in methods for nonlinear optimization problems. The algorithms used by the new Optimization Toolbox for large-scale problems include trust region methods, preconditioned conjugate gradients, projection methods for equality constraints, and reflective Newton methods for bound constraints. Again, MATLAB's sparse matrix operations form the computational core.

Trust region methods for unconstrained problems involve a local quadratic approximation to the objective function. The key components are the Hessian, a symmetric matrix of second derivatives; the gradient, a vector of first derivatives; and the trust region, an ellipsoidal subspace where the quadratic approximation is reasonably accurate. The trust region model problem is to find the vector s that minimizes

sum99_eq1.gif

subject to the constraint

sum99_eq2.gif

Here H is the Hessian, g is the gradient, D is a diagonal scaling matrix, and δ parameterizes the trust region.

For large scale problems involving thousands of variables, solving even this model problem at each step is prohibitively expensive. An additional approximation is made by restricting each minimization to the two-dimensional subspace of the trust region spanned by two vectors, s1 and s2. The vector g is the gradient s1 and the vector s2 is either the Newton direction, obtained by solving

sum99_eq3.gif
or the direction of negative curvature
sum99_eq4.gif

When the Hessian H is sparse, the direction s2 can be computed with MATLAB's sparse backslash operator or approximated by a preconditioned conjugate gradient iteration. With these choices of search directions, the method has guaranteed global convergence and fast local convergence.

Different toolbox functions take advantage of the structure of various classes of nonlinear optimization problems. There are functions for:

  • Linear equality constraints
  • Box constraints (upper and lower bounds)
  • Constrained linear least-squares
  • Nonlinear least-squares
  • Quadratic programming
  • Unconstrained nonlinear optimization
sum99_cleve_fig1.gif
sum99_cleve_fig2.gif
sum99_cleve_fig3.gif
sum99_cleve_fig4.gif

We have modified the circustent example in the toolbox to provide the example that produces the graphics in this article. Imagine constructing a circus tent to cover a square lot. The tent has five poles that will be covered with an elastic material. The cover is constrained to go over the five poles. The natural shape obtained by the tent minimizes an energy function involving its position and gradient.

The tent is approximated by specifying its height at the points in a 50-by-50 grid. These 2500 values make up the vector x. The quadratic objective function to be minimized represents the energy of the configuration. It is

sum99_eq5.gif
where H is the five-point discrete Laplacian finite difference operator (obtained from delsq and numgrid) and c is a scalar multiple of the vector of all ones. The solution must also satisfy a lower bound constraint
sum99_eq6.gif
where L is zero at most of the grid points, but has positive values at the locations of the five tent poles. The center pole is higher than the other four.
sum99_cleve_fig5.gif
sum99_cleve_fig6.gif
sum99_cleve_fig7.gif
sum99_cleve_fig8.gif

Since the objective function is quadratic and the constraints are linear, we can use the function quadprog. The first graphic shows the five tent poles and an initial configuration where all the components of x are set to the height of the center pole. The subsequent graphics show steps in the optimization. Each step reduces the energy while maintaining the constraint. After nine steps, the solution is obtained to within the accuracy of the graphics.

sum99_cleve_fig9.gif

Published 1999