vectorised matching of values with numerical tolerance

user_9319556 0
Viewed 175

I'm trying to match values in a numeric vector to a discrete set of values, but due to rounding errors they're not necessarily an exact match. Consider the following example, converting spherical coordinates to cartesian, and trying to match the resulting distance against the original r values,

phithetar <- expand.grid(theta=seq(0,pi,length=10),
                        phi = seq(0, 2*pi, length=10),
                        r = seq(1,10))

xyz <- transform(phithetar, x = r*sin(theta)*cos(phi),
                                 y = r*sin(theta)*sin(phi),
                                 z = r*cos(theta))

xyz <- transform(xyz,  newr = sqrt(x^2+y^2+z^2))

with(xyz, sum(newr != r))

This is non-zero, as the operations have results in lost precision.

Here's my slow implementation to match values within a given tolerance; the iteration over x makes it slow. An alternative would be to expand all (x,y) combinations but that'd be memory-heavy for large vectors. 

match_closest <- function(x, y, tol=1e-5){

  match_one <- function(.x) {
    difs <- abs(.x - y)
    best <- which.min(difs)
    if(difs[best] > tol) message(.x, " ain't so close to ", y[best])
    y[best]
  }
  vapply(x, match_one, 0)

}

rvalues <- unique(phithetar$r)
xyz <- transform(xyz,  matchedr = match_closest(newr, rvalues))

with(xyz, sum(matchedr != r))

I tried another idea with cut(),

match_closenough <- function(x, y, tol=1e-5){

  levs <- c(y-0.5*tol, y+0.5*tol)
  f <- cut(x, breaks=levs, include.lowest = TRUE)
  y[f]

}

but it doesn't quite work. Similarly, this rounding trick

match_closenough <- function(x, y, tol=1e-3){

  dig <- -log10(tol)
  y[match(round(x, dig), round(y, dig), nomatch = NA)]

}

kind of works and is faster, but feels too fragile.

match precision r
2 Answers

I think the cut() version can be made to work with a few amendments (though factors always seem outside my control),

match_madethecut <- function(x, y, tol=1e-5){

  breaks <- sort(c(y-0.5*tol, y+0.5*tol))
  f <- cut(x, breaks=breaks, include.lowest = TRUE)
  rep(sort(y),each=2)[as.numeric(f)]

}

Note that y is expected to contain unique values, in practice. This seems fast, as expected since cut() does the comparisons at C level.

user_93195560

If your objective is just to match values within some (very-)small tolerance to account for loss-of-precision during sequential floating-point operations (as in the example), then I think that your last match_closenough() solution is elegant and fast, and doesn't look particularly 'fragile' to me.

The only possible limitation with this approach in general seems that you're restricted to rounding by powers-of-10 (rather than arbitrarily-chosen tol values), and that you can't easily use tol values greater than 1. Both of these cases don't really apply your example, of course, but there might be situations where they might become relevant.

So: what if we insist on choosing arbitrary tol values, and want a vectorized function? One solution would be to write it in RCpp (and this would probably not be any longer than the R-only solution...).

But this can be done 'from scratch' in R: this solution uses several steps in order to avoid ever requiring a non-vectorized function:

# define sloppy_match:
sloppy_match=function(x,y,tol=.Machine$double.eps,nomatch=NA) {
    # combine x & y indices:
        v=c(x,y)
        i=c(rep(1,length(x)),rep(0,length(y)))[order(v)]
        xi=cumsum(i==1)
        yi=cumsum(i==0)
    # find y values above and below each x value
        xrank=rank(x,ties="first")
        yi_below=yi[match(xrank,xi,NA)]
        yi_above=yi_below+1
        yi_below[yi_below==0]=NA
        yi_above[yi_above>length(y)]=NA
        ysort=sort(y)   
        y_below=ysort[yi_below]
        y_above=ysort[yi_above]
    # now choose closest y value to each x:
        closest_y=y_below
        to_change=is.na(y_below)|(x-closest_y>y_above-x)&!is.na(y_above)
        closest_y[to_change]=y_above[to_change]
    # and reject any that aren't close enough
        closest_y[abs(x-closest_y)>tol]=nomatch
    return(closest_y)
}

Try it out:

vec1=sample(1:10,5000,replace=T)+rnorm(50)*1e-6
vec2=sample(1:10,5000,replace=T)

test=sloppy_match(vec1,vec2,tol=1e-3)
all(abs(test-vec1)<1e-3)
# TRUE

How does this compare in terms of speed? 

library(microbenchmark)
microbenchmark(match(vec1,vec1),
    match_closest(vec1,vec2,tol=1e-5),
    match_closenough(vec1,vec2,tol=1e-5),
    sloppy_match(vec1,vec2,tol=1e-5),
    times=10)
# Unit: microseconds
#                                       expr        min         lq        mean      median         uq        max neval
#                          match(vec1, vec1)    190.672    191.683    217.7785    199.9135    219.783    357.662    10
#     match_closest(vec1, vec2, tol = 1e-05) 186051.018 205712.022 208032.1262 209191.3980 215335.356 221682.351    10
#  match_closenough(vec1, vec2, tol = 1e-05)    511.426    513.772    531.2815    530.7850    544.416    561.935    10
#      sloppy_match(vec1, vec2, tol = 1e-05)   4110.814   4279.376   5859.7425   4301.0610   4608.116  18587.816    10

So the built-in match() is obviously fastest when the values are exactly-equal, and for close-enough matching if rounding to less-than-one powers-of-ten is acceptable (and I bet that this will be the case in 99% of circumstances), the match_closenough() is fastest.

Finally, the sloppy_match() approach to allow matching to arbitrarily-chosen tolerances is somewhat slower than these, but the benefit of vectorization is obvious: it's already almost 40x faster than the non-vectorized match_closest() approach, with only 5000 elements, and this difference becomes even-more pronounced when the number of elements increases.

So, in the end, you'll probably want to choose between the two approaches depending on the particular tolerance requirement; but both of them will be hundreds of times faster than the non-vectorized approach for large vectors.

user_134323360
Related