One of the things I enjoy about blogging is that I often learn something new. Last week I wrote about how to optimize a function that is defined in terms of an integral. While developing the program in the article, I made some mistakes that generated SAS/IML error messages. By understanding my errors, I learned something new about the default options for the QUAD subroutine in SAS/IML software.

The QUAD subroutine supports the PEAK= option, but until last week I had never used the option. (Obviously the default value works well most of the time, since I have been using the QUAD routine successfully for many years!) I discovered that the PEAK= option can be important to ensuring that a numerical integration problem converges to an accurate answer.

The default value of the PEAK= option is 1. But what does that mean? Here's what I gleaned from the documentation:

• The PEAK= option is used to scale the integrand and to determine whether the adaptive integration algorithm needs to perform additional iterations.
• The option specifies a value of x where the magnitude of the integrand f(x) is relatively large. For integrations on infinite and semi-infinite intervals, the value of the PEAK= option specifies where the function has a lot of "mass." For example, the standard normal probability density function is centered at x=0, so you can specify PEAK=0 to integrate that function on (–∞, ∞). This example illustrates why the option is named "PEAK."
• For best results, choose the PEAK= option to be in the interior of the interval of integration. Don't choose a value that is too close to the boundary. For example, on a finite domain [a, b], the midpoint of the interval is often a good choice: PEAK=(a + b)/2.
• The PEAK= value should not be a location where the integrand is zero. That is, if PEAK=v, then make sure that f(v) ≠ 0.
• In spite of its name, the value of the PEAK= option does not need to be the location of the maximum value of the integrand. For integration on a finite domain you can specify almost any value in the interval, subject to the previous rules.

The remainder of this article uses the quadratic function f(x) = 1 – x² to illustrate how you can to control the QUAD subroutine by using the PEAK= option.

Avoid zeros of the integrand

The following SAS/IML statements define the function to integrate. The integrand has a zero at the value x = 1, which is the default value for the PEAK= option. If you attempt to integrate this function on the interval [0, 2], you see the following scary-looking error message:

proc iml;
start Func(x) global(a);
   return( 1 - x##2 );         /* note that Func(1) = 0 */
finish;
 
call quad(w, "Func", {0 2});  /* default: PEAK=1 */

Convergence could not be attained over the subinterval (0 , 2 )
 
The function might have one of the following:
     1) A slowly convergent or a divergent integral.
     2) Strong oscillations.
     3) A small scale in the integrand: in this case you can change the 
        independent variable in the integrand, or
        provide the optional vector describing roughly the mean and 
        the standard deviation of the integrand.
     4) Points of discontinuities in the interior:
        in this case, you can provide a vector containing
        the points of discontinuity and try again.

The message says that the integration did not succeed, and it identifies possible problems and offers some suggestions. For this example, the integral is not divergent, oscillatory, or discontinuous, so we should look at the third item in the list. The message seems to assume that the integrand is a statistical probability distribution, which is not true for our example. For general functions, the advice could be rewritten as "... provide values for the PEAK= and SCALE= options. The SCALE= value should provide a horizontal estimate of scale; the integrand evaluated at the PEAK= value should provide a vertical estimate of scale."

For the example, –3 ≤ f(x) ≤ 1 when x is in the interval [0, 2], so the integrand has unit scale in the vertical direction. You can choose almost any value of x for the PEAK= option, but you should avoid values where the integrand is very small. The following statement uses PEAK=0.5, but you could also use PEAK= 1.5 or PEAK=0.1234.

call quad(w, "Func", {0 2}) peak=0.5;  
print w;

In a similar way, you should make sure that the integrand is defined at the PEAK= value. Avoid singular points and points of discontinuity.

Avoid the boundaries of the interval

Here's another helpful tip: choose the PEAK= option so that it is sufficiently inside the interval of integration. Avoid choosing a value that is too close to the lower or upper limit of integration. The following statement chooses a value for the PEAK= option that is just barely inside the interval [0, 2]. A warning message results:

val = 2 - 1e-7;     /* value is barely inside interval [0,2] */
call quad(w, "Func", {0 2}) peak=val;

Due to the roundoff encountered in computing FUNC
      near the point     2.00000E+00
      convergence can be attained, but only
      with an estimated error     0.00000E+00

The warning message says that a problem occurs near the point x=2. Although the routine will continue and "convergence can be attained," the routine is letting you know that it cannot perform the integration as well as it would like to. If you move the PEAK= value farther away from the limits of integration, the warning goes away and the subroutine can integrate the function more efficiently.

So there you have it. Although the default value of the PEAK= option works well for most integrals, sometimes the programmer needs to provide the QUAD routine with additional knowledge of the integrand. When you provide the PEAK= option, be sure to avoid zeros and discontinuities of the integrand, and stay away from the boundaries of the integration interval.