Chapter Contents |
Previous |
Next |
Nonlinear Optimization Examples |
All the optimization algorithms assume that f is continuous inside the feasible region. For nonlinearly constrained optimization, this is also required for points outside the feasible region. Sometimes the objective function cannot be computed for all points of the specified feasible region; for example, the function specification may contain the SQRT or LOG function, which cannot be evaluated for negative arguments. You must make sure that the function and derivatives of the starting point can be evaluated. There are two ways to prevent large steps into infeasible regions of the parameter space during the optimization process:
All the optimization techniques except the NLPNMS subroutine require continuous first-order derivatives of the objective function f. The NLPTR, NLPNRA, and NLPNRR techniques also require continuous second-order derivatives. If you do not provide the derivatives with the IML modules "grd", "hes", or "jac", they are automatically approximated by finite difference formulas. Approximating first-order derivatives by finite differences usually requires n additional calls of the function module. Approximating second-order derivatives by finite differences using only function calls can be extremely computationally expensive. Hence, if you decide to use the NLPTR, NLPNRA, or NLPNRR subroutines, you should specify at least analytical first-order derivatives. Then, approximating second-order derivatives by finite differences requires only n or 2n additional calls of the function and gradient modules. For all input and output arguments, the subroutines assume that
Note: For some difficult applications, the finite difference approximations of derivatives that are generated by default may not be precise enough to solve the optimization or least-squares problem. In such cases, you may be able to specify better derivative approximations by using a better approximation formula. You can submit your own finite difference approximations using the IML modules "grd", "hes", "jac", or "jacnlc". See Example 11.3 for an illustration.
In many applications, calculations used in the computation of f can help compute derivatives at the same point efficiently. You can save and reuse such calculations with the GLOBAL clause. As with many other optimization packages, the subroutines perform calls of the "grd", "hes", or "jac" modules only after a call of the "fun" module.
The following statements specify modules for the function, gradient, and Hessian matrix of the Rosenbrock problem:
proc iml; start F_ROSEN(x); y1 = 10. * (x[2] - x[1] * x[1]); y2 = 1. - x[1]; f = .5 * (y1 * y1 + y2 * y2); return(f); finish F_ROSEN; start G_ROSEN(x); g = j(1,2,0.); g[1] = -200.*x[1]*(x[2]-x[1]*x[1]) - (1.-x[1]); g[2] = 100.*(x[2]-x[1]*x[1]); return(g); finish G_ROSEN; start H_ROSEN(x); h = j(2,2,0.); h[1,1] = -200.*(x[2] - 3.*x[1]*x[1]) + 1.; h[2,2] = 100.; h[1,2] = -200. * x[1]; h[2,1] = h[1,2]; return(h); finish H_ROSEN;
The following statements specify a module for the Rosenbrock function when considered as a least-squares problem. They also specify the Jacobian matrix of the least-squares functions.
proc iml; start F_ROSEN(x); y = j(1,2,0.); y[1] = 10. * (x[2] - x[1] * x[1]); y[2] = 1. - x[1]; return(y); finish F_ROSEN; start J_ROSEN(x); jac = j(2,2,0.); jac[1,1] = -20. * x[1]; jac[1,2] = 10.; jac[2,1] = -1.; jac[2,2] = 0.; return(jac); finish J_ROSEN;
Example 22 of Mor, Garbow, and Hillstrom (1981) illustrates the sparse Hessian module input. The objective function, which is the Extended Powell's Singular Function, for n=40 is a least-squares problem:
start f_nlp22(x); n=ncol(x); f = 0.; do i=1 to n-3 by 4; f1 = x[i] + 10. * x[i+1]; r2 = x[i+2] - x[i+3]; f2 = 5. * r2; r3 = x[i+1] - 2. * x[i+2]; f3 = r3 * r3; r4 = x[i] - x[i+3]; f4 = 10. * r4 * r4; f = f + f1 * f1 + r2 * f2 + f3 * f3 + r4 * r4 * f4; end; f = .5 * f; return(f); finish f_nlp22; start g_nlp22(x); n=ncol(x); g = j(1,n,0.); do i=1 to n-3 by 4; f1 = x[i] + 10. * x[i+1]; f2 = 5. * (x[i+2] - x[i+3]); r3 = x[i+1] - 2. * x[i+2]; f3 = r3 * r3; r4 = x[i] - x[i+3]; f4 = 10. * r4 * r4; g[i] = f1 + 2. * r4 * f4; g[i+1] = 10. * f1 + 2. * r3 * f3; g[i+2] = f2 - 4. * r3 * f3; g[i+3] = -f2 - 2. * r4 * f4; end; return(g); finish g_nlp22;
You can specify the sparse Hessian with the following module:
start hs_nlp22(x); n=ncol(x); nnz = 8 * (n / 4); h = j(nnz,3,0.); j = 0; do i=1 to n-3 by 4; f1 = x[i] + 10. * x[i+1]; f2 = 5. * (x[i+2] - x[i+3]); r3 = x[i+1] - 2. * x[i+2]; f3 = r3 * r3; r4 = x[i] - x[i+3]; f4 = 10. * r4 * r4; j= j + 1; h[j,1] = i; h[j,2] = i; h[j,3] = 1. + 4. * f4; h[j,3] = h[j,3] + 2. * f4; j= j+1; h[j,1] = i; h[j,2] = i+1; h[j,3] = 10.; j= j+1; h[j,1] = i; h[j,2] = i+3; h[j,3] = -4. * f4; h[j,3] = h[j,3] - 2. * f4; j= j+1; h[j,1] = i+1; h[j,2] = i+1; h[j,3] = 100. + 4. * f3; h[j,3] = h[j,3] + 2. * f3; j= j+1; h[j,1] = i+1; h[j,2] = i+2; h[j,3] = -8. * f3; h[j,3] = h[j,3] - 4. * f3; j= j+1; h[j,1] = i+2; h[j,2] = i+2; h[j,3] = 5. + 16. * f3; h[j,3] = h[j,3] + 8. * f3; j= j+1; h[j,1] = i+2; h[j,2] = i+3; h[j,3] = -5.; j= j+1; h[j,1] = i+3; h[j,2] = i+3; h[j,3] = 5. + 4. * f4; h[j,3] = h[j,3] + 2. * f4; end; return(h); finish hs_nlp22; n = 40; x = j(1,n,0.); do i=1 to n-3 by 4; x[i] = 3.; x[i+1] = -1.; x[i+3] = 1.; end; opt = j(1,11,.); opt[2]= 3; opt[9]= 8 * (n / 4); call nlpnra(xr,rc,"f_nlp22",x,opt) grd="g_nlp22" hes="hs_nlp22";
Note: If the sparse form of Hessian defines a diagonal matrix (that is, i = j in all nn rows), the NLPNRA algorithm stores and processes a diagonal matrix G. If you do not specify any general linear constraints, the NLPNRA subroutine uses only order n memory.
Chapter Contents |
Previous |
Next |
Top |
Copyright © 1999 by SAS Institute Inc., Cary, NC, USA. All rights reserved.