Re: Roundoff error in SMOOTH

[Boyd Blackwell <Boyd.Blackwell(at)anu.edu.au>]

THe IDL smooth function scales well with the number of points in the average
(i.e. smooth(x,10) and smooth(x,1000) take the same time.  It is very likely
that it uses an incremental algorithm, whereby result(i) = result(i-1) + x(i+n)
- x(i-1)
where n is the number of points in the average.
This is inherently subject to roundoff error, albeit relatively small.

I can be sure it doe not use FFTs, as it is way too fast for that.



Wayne Landsman wrote:
> 
> I was recently surprised to discover that applying the SMOOTH function to a
> non-negative array could  yield an array with negative numbers.  I give an
> example below.  This problem is evidently due to some sort of roundoff error,
> since it does not occur when using double precision.   But it is not obvious to
> me how averaging 9 non-negative numbers (for a 3x3 box smooth) could yield a
> negative number, even allowing for roundoff error.
> 
> Although not obvious from my simple example, this has nothing to do with edge
> effects -- I originally found the problem when 3x3 smoothing a 1024 x 1024
> array.
> 
> --Wayne Landsman                             landsman@mpb.gsfc.nasa.gov
> 
> d = fltarr(11,11)                 ;Initialize array to zero
> d[ [0,2,0,2],[2,3,4,5] ] = 0.48   ;Set a few values to a positive number
> print,where(d lt 0)               ;Make sure there are really no negative values
> ;          -1
> 
> print,d[1:6,1:6]                            ;Print a block of numbers in array
> ;      0.00000      0.00000      0.00000      0.00000      0.00000      0.00000
> ;      0.00000      0.00000      0.00000      0.00000      0.00000      0.00000
> ;      0.00000     0.480000      0.00000      0.00000      0.00000      0.00000
> ;      0.00000      0.00000      0.00000      0.00000      0.00000      0.00000
> ;     0.480000      0.00000      0.00000      0.00000      0.00000      0.00000
> ;      0.00000      0.00000      0.00000      0.00000      0.00000      0.00000
> 
> dd = smooth(d,3)       ;Smooth the array
> print,dd[1:6,1:6]      ;After smooothing there are negative values!
> ;      0.00000      0.00000      0.00000      0.00000      0.00000      0.00000
> ;    0.0533333    0.0533333    0.0533333      0.00000      0.00000      0.00000
> ;     0.106667    0.0533333    0.0533333      0.00000      0.00000      0.00000
> ;     0.160000     0.106667    0.0533333 -6.62274e-09 -6.62274e-09 -6.62274e-09
> ;     0.106667    0.0533333      0.00000      0.00000      0.00000      0.00000
> ;    0.0533333    0.0533333      0.00000      0.00000      0.00000      0.00000
> 
> dd = smooth(double(d),3)   ;Negative values do not occur in double precision
> print,float(dd[1:6,1:6])
> ;      0.00000      0.00000      0.00000      0.00000      0.00000      0.00000
> ;    0.0533333    0.0533333    0.0533333      0.00000      0.00000      0.00000
> ;     0.106667    0.0533333    0.0533333      0.00000      0.00000      0.00000
> ;     0.160000     0.106667    0.0533333      0.00000      0.00000      0.00000
> ;     0.106667    0.0533333      0.00000      0.00000      0.00000      0.00000
> ;    0.0533333    0.0533333      0.00000      0.00000      0.00000      0.00000