eddelbuettel / mkl4deb

Adding the Intel MKL to a Debian / Ubuntu system via one simple script
212 stars 37 forks source link

MKL inconsistent results #13

Closed CYGUBICKO closed 3 years ago

CYGUBICKO commented 3 years ago

I have R package pcoxtime (my implementation but already in CRAN) which I am using to fit survival model using proximal gradient descent. I have noticed that whenever I run the code on a system which has MKL as the default, I get different (seemingly wrong) results. This mostly happens when the dataset is large. For small dataset (few observations), the results are comparable (generally similar) with other packages (what we expect). Here's part of my sessionInfo() in this case.

R version 3.6.3 (2020-02-29)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.5 LTS

Matrix products: default
BLAS/LAPACK: /opt/intel/compilers_and_libraries_2018.2.199/linux/mkl/lib/intel64_lin/libmkl_rt.so

However, default BLAS/LAPACK libraries give me expected results which are consistent. No matter the size of the data.

One thing to mention, pcoxtime has most of it's matrix operations implemented using Rcpp and RcppArmadillo.

Here is the example, but first install pcoxtime and also the problem only occurs if you MKL as the default.

Small n: print(mod_coxph) and print(mod_pcox) are close enough, which is what we expect.

# install.packages("pcoxtime")
library(survival)
library(pcoxtime)

set.seed(911)
sims <- simtdc(10)
df <- sims$data
betas <- sims$betas
print(betas)
df <- df[, c("Start", "Stop", "Event", names(betas))]
## coxph
mod_coxph <- coxph(Surv(Start, Stop, Event) ~ ., data = df, method = "breslow")
print(mod_coxph)
## pcoxtime
mod_pcox <- pcoxtime(Surv(Start, Stop, Event) ~ ., data = df, alpha = 1, lambda = 0)
print(mod_pcox)

Relatively large n: print(mod_coxph) and print(mod_pcox) gives inconsistent results

sims <- simtdc(100)
df <- sims$data
betas <- sims$betas
print(betas)
df <- df[, c("Start", "Stop", "Event", names(betas))]
## coxph
mod_coxph <- coxph(Surv(Start, Stop, Event) ~ ., data = df, method = "breslow")
print(mod_coxph)
## pcoxtime
mod_pcox <- pcoxtime(Surv(Start, Stop, Event) ~ ., data = df, alpha = 1, lambda = 0)
print(mod_pcox)

Any help will be appreciated.

Thanks,

eddelbuettel commented 3 years ago

I don't represent Intel or the MKL team; I simply provided a helper script for the platforms I work on / with. So I think we should close this here.

If you have an identifiable issue you could possibly test, and then warn, from your package's configure script. Otherwise there is zero we can do: as R package authors we have no influence over what BLAS / LAPACK libraries we will find with R.