4969 Rate this article:
No rating

Beware the behavior of SMOOTH

Mark Piper

Suppose I have an array that has an outlier. A really big outlier:

IDL> a = [1.0, 1.0, 2.0, 3.0, 4.0, 1.0d18, 4.0, 3.0, 2.0, 1.0, 1.0]
I'd like to smooth this double-precision array with a running mean (or boxcar or tophat, depending on where you learned this technique) filter, such as provided by the IDL SMOOTH function. For efficiency, SMOOTH divides the kernel width into a running total of the differences between neighboring values. Here's the result of applying SMOOTH to A with a filter width of 3:
IDL> print, smooth(a, 3)
       1.0000000       1.3333333       2.0000000       3.0000000  3.3333333e+017
  3.3333333e+017  3.3333333e+017      0.00000000      -1.0000000      -1.6666667
Whoa. The input array is symmetric, so why isn't the output? Also—and this is worrisome—the input array is composed of positive numbers, so how can the mean of any subset of these numbers be negative? The answer lies in the way floating point numbers are represented on computers. To see why, we can use the information returned from the MACHAR function:
IDL> m = machar(/double)
IDL> help, m
** Structure DMACHAR, 13 tags, length=72, data length=68:
IBETA           LONG                 2
IT              LONG                53
IRND            LONG                 5
NGRD            LONG                 0
MACHEP          LONG               -52
NEGEP           LONG               -53
IEXP            LONG                11
MINEXP          LONG             -1022
MAXEXP          LONG              1024
EPS             DOUBLE      2.2204460e-016
EPSNEG          DOUBLE      1.1102230e-016
XMIN            DOUBLE      2.2250739e-308
XMAX            DOUBLE      1.7976931e+308
Using the IBETA field, which gives the base used to construct numbers (2, no surprise), and the IT field, which gives the number of base-2 digits used in the mantissa of a number, the maximum resolvable distance between two double precision numbers must be given by:
IDL> mrd = double(m.ibeta)^m.it
IDL> print, mrd
So, when SMOOTH tries to difference two numbers whose distance is greater than MRD, bad things can happen because of loss of precision. Let's apply this information to the example above. Define values slightly above and below the threshold set by MRD:
IDL> below = double(m.ibeta)^(m.it-1)
IDL> above = double(m.ibeta)^(m.it+1)
IDL> print, below, above
  4.5035996e+015  1.8014399e+016
and substitute them into the array used above:
IDL> b = [1.0, 1.0, 2.0, 3.0, 4.0, below, 4.0, 3.0, 2.0, 1.0, 1.0]
IDL> c = [1.0, 1.0, 2.0, 3.0, 4.0, above, 4.0, 3.0, 2.0, 1.0, 1.0]
Now apply SMOOTH to these arrays and evaluate the results:
IDL> print, smooth(b, 3)
       1.0000000       1.3333333       2.0000000       3.0000000  1.5011999e+015
  1.5011999e+015  1.5011999e+015       3.0000000       2.0000000       1.3333333
IDL> print, smooth(c, 3)
       1.0000000       1.3333333       2.0000000       3.0000000  6.0047995e+015
  6.0047995e+015  6.0047995e+015       3.3333333       2.3333333       1.6666667
Note that SMOOTH works well when applied to B, but not to C—the results aren't symmetric.

A better discussion of this behavior is given in the section Note on Smoothing Over Large Data Ranges in the IDL Help page for SMOOTH, along with a workaround for this situation.

Update (2013-07-01): I neglected to mention that the seed for this post, as well as the note in the Help, came from discussions a few years ago with Carmen Lucas at DRDC Atlantic. Thanks, Carmen, for pointing out this unexpected behavior to me.


A Technical Dive into Creating Accurate Digital Terrain Models


Live Webinar | July 25 Creating accurate digital terrain models (DTMs) is essential for modeling areas of... more »

The Science, and Art, of Hyperspectral Image Analysis


Harris Geospatial Solutions – June 7, 2018 Automation is being touted as the wave of the future for... more »


Sign up to receive the latest news, events, technologies and special offers.