1
$\begingroup$

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.

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

The various Fortran standards allow a lot of compiler dependent behaviour in terms of function binary interfaces when being called with "complicated" data types such as Fortran90 style arrays and complex numbers. This means calling code compiled with one compiler from another one is not guaranteed to do what you expect, and can lead to grabbing the wrong bit of memory for a variable, thus giving spurious results.

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.

$\endgroup$
  • $\begingroup$ Thanks - that's also what I was thinking and having a second supporting opinion definitely helps ;) $\endgroup$ – xeroqu Aug 25 '17 at 12:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.