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.
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.
1 Answer
Reset to default 0Short 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./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/content/www/us/en/developer/articles/technical/introduction-to-the-conditional-numerical-reproducibility-cnr.html)