R matrix calculations less precise with Intel oneMKL

Question

I have installed Intel oneAPI Math Kernel Library (oneMKL) on Windows 10 and then used the trick described here to slip those DLL's to R without having to recompile R. The 7-fold speedup of matrix operations is amazing. However, I observed some numerical imprecisions that didn't occur before.

In particular, I have a code that uses matrix multiplications to calculate a matrix of squared euclidean distances. There should be obviously zeros on the diagonal, but there are sometimes even negative values (which causes problems when you do sqrt). The largest imprecision is -7.275958e-12. Isn't it quite large, given the precision of double is 1e-15?

The code below demonstrates it on part of my data:

# locations
x1 <- structure(c(-731.312436904721, -762.671574325304, -696.47364436844, 
-703.029671223004, -672.797495061876, -697.913648506778, -692.427540613256, 
-666.452378882586, -648.220100118672, -604.945020205791, -946.343545843331, 
-958.781357545289, -973.468353412248, -978.218709475081, -973.701316142042, 
-984.499258109887, -988.016827462228, -994.1159260048, -990.748066075698, 
-1001.40046471906), dim = c(10L, 2L), dimnames = list(NULL, c("x_km", 
"y_km")))

n1 <- nrow(x1)

# my matrix trick for fast computation of matrix of squared euclidean distances
x1 <- x1 / rep(c(10,10), each = n1) # scaling
x1_2 <- apply(x1^2, 1, sum)
sq_dist <- cbind(x1, x1_2, 1) %*% rbind(t(-2*x1), 1, x1_2)


range(diag(sq_dist))
# [1] -0.000000000003637979  0.000000000003637979

So, I will probably have to do something like sq_dist[sq_dist < 0] <- 0, but I also wonder whether this isn't really a bug, because the imprecision of ~7e-12 seems quite large given the precision of double around 1e-15.

Kyler · Accepted Answer · 2025-03-14 08:34:10Z

1

Short answer: This is probably nothing major to worry about and is expected behavior using those specific optimizations. Rounding or driving very small numbers to zero is your best bet.

Long answer: This is a common floating point error (ref: https://floating-point-gui.de/). You are likely getting similar results normally but R is kind enough to round them for you in a standard environment. If you want to use scipen you can disable printing scientific notation (ref: https://cran.r-project.org/web/packages/Tplyr/vignettes/options.html#:~:text=tplyr.scipen,formatted%20in%20presentation.). Also, reading the developer guide for that oneMKL, you may also be able to at least make it produce more consistent results using `MKL_CBWR = AVX2,STRICT` (ref: https://www.intel.com/content/www/us/en/developer/articles/technical/introduction-to-the-conditional-numerical-reproducibility-cnr.html)

edited Mar 14 at 8:34

answered Mar 14 at 8:29

Kyler

113 bronze badges

Sign up to request clarification or add additional context in comments.

2 Comments

Tomas Mar 14 at 11:48

Thanks for your response! 1) Well, the "common floating error" you reference is 1e-17, this is 1e-12, that's 10000x larger.. 2) I don't think "R is kind enough to round this for me". How would R know that there should be zeros on the diagonal? It can't.

Tomas Mar 14 at 11:56

and thanks for the info on MKL_CBWR setting! I tried Sys.setenv(MKL_CBWR = "AVX2,STRICT"), and also setting it before calling R, but the result is still the same.

Collectives™ on Stack Overflow

R matrix calculations less precise with Intel oneMKL

1 Answer 1

2 Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

2 Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related