# The difference between mkl_intel_lp64 vs mkl_gf_lp64 in a numerical reproducibility issue with Intel MKL

It has been discussed that Intel MKL can exhibit irreproducible behavior under certain conditions. In fact, this is a known thing and described by Intel as Conditional Numerical Reproducibility.

A colleague has recently encountered this issue while trying to find the eigenvalues of a matrix and he realized that the function eigen in R is not deterministic on the server he was using. In particular, eigenvalues can be slightly different (<1e-12) and the signs of the eigenvectors are non-deterministic.

I tested his hypothesis by considering this matrix

 corr <- structure(c(1, 0.250050163743823, 0.23435347961023, 0.482039807584937,
0.260244618444847, 0.245304368749426, 0.486808023309575, 0.479348415857738,
0.250050163743823, 1, 0.24648658148082, 0.246713558484002, 0.249547953702824,
0.249424786521267, 0.247921981042718, 0.493767069248003, 0.23435347961023,
0.24648658148082, 1, 0.0140300592972217, 0.229687680789224, 0.251843633087952,
0.0102624595100619, 0.0160021143836412, 0.482039807584937, 0.246713558484002,
0.0140300592972217, 1, 0.47661703967983, 0.00364166259778339,
0.257242621615393, 0.507588614471663, 0.260244618444847, 0.249547953702824,
0.229687680789224, 0.47661703967983, 1, 0.243407621930236, 0.482342591096063,
0.473425275160343, 0.245304368749426, 0.249424786521267, 0.251843633087952,
0.00364166259778339, 0.243407621930236, 1, 0.00266373582273809,
0.00415353194451873, 0.486808023309575, 0.247921981042718, 0.0102624595100619,
0.257242621615393, 0.482342591096063, 0.00266373582273809, 1,
0.505031653312059, 0.479348415857738, 0.493767069248003, 0.0160021143836412,
0.507588614471663, 0.473425275160343, 0.00415353194451873, 0.505031653312059,
1), .Dim = structure(c(8L, 8L), .Names = c("rt", "rt")), .Dimnames = structure(list(
rt = c("A", "B", "C", "D", "E", "F", "G", "H"
), rt = c("A", "B", "C", "D", "E", "F", "G",
"H")), .Names = c("rt", "rt")))


after computing its eigenvectors and eigenvalues I saw that they are indeed not the same:

> percDet <- 100 * mean(sapply(1:1000, function(...) {
+     identical(eigen(corr, symmetric = TRUE),
+               eigen(corr, symmetric = TRUE))}))
> message("% eigen determinism: ", percDet, "%")
% eigen determinism: 75.7%

> percDet <- 100 * mean(sapply(1:1000, function(...) {
+     identical(eigen(corr, symmetric = TRUE)$values, + eigen(corr, symmetric = TRUE)$values)}))
> message("% eigenvalues determinism: ", percDet, "%")
% eigenvalues determinism: 76.3%


After some investigation I've found out that the issue was due to the Intel MKL configuration we used for R. Specifically, in our Ansible role we were configuring R as follows:

./configure --enable-R-shlib --with-blas=\"-Wl,--no-as-needed -L${MKLROOT}/lib/intel64 -L{{ mkl_install_dir }}/compiler/lib/intel64 -lmkl_gf_lp64 -lmkl_core -lmkl_intel_thread -liomp5 -lpthread -lm\"  Changing "-lmkl_gf_lp64" part to "-lmkl_intel_lp64" as below ./configure --enable-R-shlib --with-blas="-Wl,--no-as-needed -L${MKLROOT}/lib/intel64 -L{{ mkl_install_dir }}/compiler/lib/intel64 -lmkl_intel_lp64 -lmkl_core -lmkl_intel_thread -liomp5 -lpthread\"


solved the issue:

> percDet <- 100 * mean(sapply(1:1000, function(...) {
+ identical(eigen(corr, symmetric = TRUE),
+ eigen(corr, symmetric = TRUE))}))
> message("% eigen determinism: ", percDet, "%")
% eigen determinism: 100%

> percDet <- 100 * mean(sapply(1:1000, function(...) {
+ identical(eigen(corr, symmetric = TRUE)$values, + eigen(corr, symmetric = TRUE)$values)}))
> message("% eigenvalues determinism: ", percDet, "%")
% eigenvalues determinism: 100%


While I have found a related question "Intel MKL - Difference between mkl_intel_lp64 and mkl_gf_lp64" helpful, it doesn't explain why this leads to a reproducibility issue when multiple threads are used.

So, I'm wondering whether there is a good explanation to this behavior and why changing "mkl_gf_lp64" flag to "mkl_intel_lp64" resolves the issue.

• I don't believe that you've turned on the CNR feature in either case, so it's possibly coincidence that this change worked. – Brian Borchers Sep 4 '17 at 20:11

You're using the intel compiler suite to build R, so you need to to link to the intel build of the fortran library for your BLAS to behave correctly. This is what the -lmkl_intel_lp64, embedded inside an instruction to the linker, does (as opposed to your previous version, which was linking against a gfortran version, `libmkl_gf_lp64"). Since threading makes it more likely that the wrong memory you grabbed was different between different runs, this showed up as non-deterministic behaviour, but there's no prior reason to thing that a serial run linking against the wrong library was generating the right result either.