APPLICATIONS OF SAS IML: SPH 5421 notes.018
FINDING MAXIMA AND MINIMA OF FUNCTIONS USING IML
Newton's method gives a way of finding a solution to a nonlinear vector
equation of the form
G(W) = 0,
where W is a vector in R^n, and 0 is the 0 vector in R^m.
Assume H(W) is a function from a subset S of R^n to R^1. That is, it
is a real-valued function of an n-vector. It is often desirable to find the
minimum value of such a function. This means, find the value W0 of
W such that H(W0) is smaller than H(W) for any other W in the subset S.
Here is a simple example. Assume W is a two-dimensional vector,
W = (u, v)'. Assume H(W) = H(u, v)' = u^2 + 5 * v^2 + 12. Then H assumes its
smallest value at the point u = v = 0. You don't need any complex computations
to find this.
In general it is harder to find where a given function might take on its
minimum value. If it is a smooth (differentiable) function which takes on its
minimum on the interior of the subset S where it is defined, then all the
partial derivatives at that point will be 0. The problem of finding where the
minimum occurs is thus reduced to finding the points where all the partial
derivatives are zero.
Newton's method thus provides a way of locating minima (or maxima) of a
differentiable function. One simply takes the partial derivatives of the
function with respect to the parameters, and sets them equal to zero. If the
original function is H(W) = H(u, v), the equations are
dH(u, v)/du = 0
dH(u, v)/dv = 0
This can be abbreviated to:
dH(u, v) = 0,
where 0 is actually the vector (0, 0)'.
To apply Newton's method to this system of nonlinear equations, you
need to compute the SECOND derivatives of H. The resulting matrix is
| d2H/du2 d2H/dudv |
d2H = | |
| d2H/dudv d2H/dv2 |
It is now straightforward to apply Newton's method. Assume w0 = (u0, v0)' is
an initial guess at a solution. Then define
w1 = w0 - inv(d2H(w0)) * dH(w0)
Continue iterations in this way until criteria for convergence are met.
Convergence criteria might be stated in terms of the distance between
successive iterations:
dist(w(i+1), wi) = max{u(i+1) - yi, v(i+1), vi}
If the distance is smaller than some pre-defined level, (for example,
1e-6), the iterations are stopped.
The following is a SAS program which implements this algorithm
using PROC IML, and printout from the program:
========================================================================
/* Program to compute minima or maxima of a nonlinear real-valued */;
/* function - based on Newton-Raphson methods - uses IML */;
footnote "program: /home/walleye/john-c/5421/imlmin.sas &sysdate &systime" ;
options linesize = 80 ;
proc iml ;
beta = { .5, .5} ; /* Initialize the param. vector */;
/* -------------------------------------------------------------------*/;
/* Module which computes the function and its derivatives ... */;
/* -------------------------------------------------------------------*/;
start funct(beta, h, dh, d2h) ;
x = beta[1] ;
y = beta[2] ;
h = exp(x*x + 3) + exp(y*y - 5) ; /* function defn */
dhdx = 2*x*exp(x*x + 3) ; /* 1st derivs */
dhdy = 2*y*exp(y*y - 5) ;
d2hdx2 = 2*exp(x*x + 3) + 4*x*x*exp(x*x + 3) ; /* 2nd derivs */
d2hdxdy = 0 ;
d2hdy2 = 2*exp(y*y - 5) + 4*y*y*exp(y*y - 5) ;
dh = j(2, 1, 0) ; /* Initialize 1st deriv vector */
d2h = j(2, 2, 0) ; /* Initialize 2nd deriv matrix */
dh[1] = dhdx ;
dh[2] = dhdy ;
d2h[1, 1] = d2hdx2 ;
d2h[1, 2] = d2hdxdy ;
d2h[2, 1] = d2hdxdy ;
d2h[2, 2] = d2hdy2 ;
finish funct ;
/* -------------------------------------------------------------------*/;
eps = 1e-8 ;
diff = 1 ;
/* The following do loop stops when the increments in beta are */;
/* sufficiently small, or when number of iterations reaches 30. */;
do iter = 1 to 30 while(diff > eps) ;
run funct(beta, h, dh, d2h) ;
beta = beta - inv(d2h) * dh ; /* The key Newton step ... */;
diff = max(abs(inv(d2h) * dh)) ;
print iter h dh d2h diff beta ;
end ;
========================================================================
17:27 Tuesday, February 15, 2000
ITER H DH D2H DIFF BETA
1 25.798992 25.79034 77.37102 0 0.3333333 0.1666667
0.0086517 0 0.0259551 0.1666667
ITER H DH D2H DIFF BETA
2 20.658218 6.8837633 43.597167 0 0.1578947 0.0087719
0.0023092 0 0.0146252 0.0087719
ITER H DH D2H DIFF BETA
3 20.093821 0.352405 40.180348 0 0.0087706 1.3497E-6
0.0001182 0 0.013479 1.3497E-6
ITER H DH D2H DIFF BETA
4 20.092275 0.0000542 40.171074 0 1.3497E-6 4.918E-18
1.8189E-8 0 0.0134759 4.918E-18
ITER H DH D2H DIFF BETA
5 20.092275 1.976E-16 40.171074 0 4.918E-18 0
6.627E-20 0 0.0134759 0
program: /home/walleye/john-c/5421/immin.sas 15FEB00 17:27
========================================================================
NOTES ON THE PRECEDING PROGRAM / PRINTOUT
1. You can tell by inspection that the function has a minimum value
at (0, 0)'. So it is no surprise that that is what the program
finds.
2. Note that the value of the function, H, decreases from one iteration
to the next. The algorithm appears to be finding a relative minimum.
3. Note that both components of the 1st derivative vector, DH, get
smaller as you go from one iteration to the next. The 2 x 2 matrix
of 2nd derivatives, D2H, converges to a nonzero constant. Note
that it has positive determinant, as would be expected for a
relative minimum.
========================================================================
PROBLEM 17
1. For the function G(x, y) = x^4 + y^4 - 5*x^2 -7*y^2 + 12, write
a SAS program that finds where G(x, y) takes on its minimum value.
2. Use R or SAS to draw a graph or contour plot of the function G(x, y).
========================================================================
/home/walleye/john-c/5421/notes.018 Last update: November 9, 2012