[R] MASS::mvrnorm() on MKL may produce different numbers even when the seed is the same?
Shu Fai Cheung
@hu|@|@cheung @end|ng |rom gm@||@com
Thu Aug 17 16:33:17 CEST 2023
Sorry for the confusion caused by that sentence. I did
not expect using the same seed generates the same
sequence of numbers unconditionally. What puzzled me
is not having the same sequence when running the
same code repeatedly in the same computer in the same
session of R. The following four sets of results were generated
from a script with four copies of the two lines, ran in the
same session:
> set.seed(1234)
> MASS::mvrnorm(n = 5, mu = mu0, Sigma = sigma0)[1, 1:5]
[1] 0.4851605 0.5704446 1.6873036 2.7645014 -0.2020908
> set.seed(1234)
> MASS::mvrnorm(n = 5, mu = mu0, Sigma = sigma0)[1, 1:5]
[1] 0.4851605 0.5704446 1.6872818 2.7644787 -0.2023346
> set.seed(1234)
> MASS::mvrnorm(n = 5, mu = mu0, Sigma = sigma0)[1, 1:5]
[1] 0.4851605 0.5704446 1.6872818 2.7644787 -0.2023346
> set.seed(1234)
> MASS::mvrnorm(n = 5, mu = mu0, Sigma = sigma0)[1, 1:5]
[1] 0.4851605 0.5704446 1.6873036 2.7645014 -0.2020908
I could not figure out the pattern. I could not even
reproduce the order above. Sometimes, the order can be
different. Sometimes, the four sets of results can be identical.
I did a few more tests, on another computer (Windows 10),
with several copies of R installed. For R 4.3.1 and
R-devel (2023-08-12 r84939 ucrt), for which I did not touch
anything about BLAS and LAPACK, they consistently gave
the following results:
> set.seed(1234)
> MASS::mvrnorm(n = 5, mu = mu0, Sigma = sigma0)[1, 1:5]
[1] 0.4851605 0.5704446 1.6872818 2.7644787 -0.2023346
Therefore, it is possible that in my test with MKL, somehow
BLAS and LAPACK came with R were sometimes used, even
though the sessionInfo() output shows that MKL is used.
I am not sure whether it is because I did anything wrong in
setting up R to use MKL.
More on the machine on which I found that issue. It was a virtual
machine built using VirtualBox. The host system is Windows 10
(though I guess it should not matter). The OS is Ubuntu 22.04.3 LTS.
I installed MKL (2023.2.0) following these steps:
https://www.intel.com/content/www/us/en/developer/tools/oneapi/onemkl-download.html?operatingsystem=linux&distributions=aptpackagemanager
I compiled R-devel (2023-08-15 r84957) and then link MKL to it
following these steps (with some modifications, as the description
seems to be outdated):
https://www.intel.com/content/www/us/en/developer/articles/technical/quick-linking-intel-mkl-blas-lapack-to-r.html
Ideally, I should compile R with MKL, instead of using the quick
linking method. However, I am occupied with something else and
could not find an updated guide on how to do this quickly. The
guides I located were several years old.
Regards,
Shu Fai
On Thu, Aug 17, 2023 at 8:57 PM Ben Bolker <bbolker using gmail.com> wrote:
>
> > However, should the numbers
> > generated identical if the same seed is used?
>
> I don't see how using the same seed can overcome floating-point
> differences across platforms (compilers etc.) stemming from differences
> in an eigen() computation (based on arcane details like use of
> registers, compiler optimizations, etc.). set.seed() guarantees that
> the normal deviates generated by rnorm() will be identical, but those
> random numbers are then multiplied by the eigenvectors derived from
> eigen(), at which point the differences will crop up.
>
> This has been discussed on Twitter:
> https://twitter.com/wviechtb/status/1230078883317387264
>
> Wolfgang Viechtbauer, 2020-02-18:
>
> Another interesting reproducibility issue that came up. MASS::mvrnorm()
> can give different values despite setting the same seed. The problem:
> Results of eigen() (which is used by mvrnorm) can be machine dependent
> (help(eigen) does warn about this).
>
> Interestingly, mvtnorm::rmvnorm() with method="eigen" gives the same
> values across all machines I tested this on (and method="svd" give the
> same values). With method="chol" I get different values, but again
> consistent across machines.
>
> Ah, mvtnorm::rmvnorm() applies the results from eigen() in a different
> way that appears to be less (not?) affected by the indeterminacy in the
> eigenvectors. Quite clever.
>
>
> On 2023-08-16 11:10 p.m., Shu Fai Cheung wrote:
> > Hi All,
> >
> > When addressing an error in one of my packages in the CRAN check at CRAN,
> > I found something strange which, I believe, is unrelated to my package.
> > I am not sure whether it is a bug or a known issue. Therefore, I would
> like
> > to have advice from experts here.
> >
> > The error at CRAN check occurred in a test, and only on R-devel with
> MKL. No
> > issues on all other platforms. The test compares two sets of random
> numbers,
> > supposed to be identical because generated from the same seed. However,
> when
> > I tried that in R-devel (2023-08-15 r84957) on Ubuntu 22.04.3 LTS,
> linked to
> > MKL (2023.2.0), I found that this is not the case in one situation.
> >
> > I don't know why but somehow only one particular set of means and
> > covariance matrix, among many others in the code, led to that problem.
> > Please find
> > below the code to reproduce the means and the covariance matrix (pardon
> me
> > for the long code):
> >
> > mu0 <- structure(c(0.52252416853188655, 0.39883382931927774,
> 1.6296642535174521,
> > 2.9045763671121816, -0.19816874840500939, 0.42610841566522556,
> > 0.30155498531316366, 1.0339601619394503, 3.4125587827873192,
> > 13.125481598475405, 19.275480386183503, 658.94225353462195,
> 1.0997193726987393,
> > 9.9980286642877214, 6.4947188998602092, -12.952617780813679,
> > 8.4991882784024799, 106.82209520278377, 0.19623712449219838), class =
> > c("numeric"))
> >
> > sigma0 <- structure(c(0.0047010182215945009, 1.4688424622163808e-17,
> > -2.5269697696885822e-17,
> > -1.2451516929358655e-18, 3.6553475725253408e-19, 3.2092363914356227e-19,
> > -2.6341938382508503e-18, 8.6439582878287556e-19, -1.295240602054123e-17,
> > 1.1154057497258702e-18, -8.2407621704807367e-33, -1.0891686096304908e-15,
> > -2.6800671450262819e-18, -0.00092251429799113333,
> -1.4833018148344697e-16,
> > 2.2888892242132203e-16, -3.5128615476876035e-18, 2.8004770623734002e-33,
> > 2.0818099301348832e-34, 1.294907062174661e-17, 0.012788608283099183,
> > 1.5603789410838518e-15, 1.0425182851251724e-16, 2.758318745465066e-17,
> > -9.7047963858148572e-18, -7.4685438142667792e-17,
> -1.9426723638193857e-17,
> > -7.8775307262071738e-17, -4.0773523140359949e-17, 4.1212975037936416e-18,
> > -1.199934828824835e-14, 2.301740948367742e-17, -2.5410883836579325e-18,
> > -0.129172198695584, -1.9922957470809647e-14, 1.7077690739402761e-16,
> > 1.663702957392817e-30, 5.5362608501514495e-32, -2.5575000656645995e-17,
> > 7.3879960338026333e-16, 0.052246980157830372, -0.007670030643844231,
> > -1.5468402204507243e-17, 1.0518757786432683e-17, 2.9992131812421253e-16,
> > 1.6588830885081527e-17, -8.3276045185697992e-16, -3.3713008849128695e-15,
> > 1.5823980389961169e-17, 2.2243227896948194e-15, -3.2071920299461956e-17,
> > 5.0187645877462975e-18, -7.4622951185527265e-15, -0.44701112744184102,
> > -1.854066605407503e-16, -5.2511356793649759e-30, -6.7993385117653912e-32,
> > -1.5343740968067595e-18, -3.6785921616583949e-17, -0.0076700306438433133,
> > 0.019231143599164325, 5.7818784939170702e-18, 1.7819897158350214e-17,
> > -2.3927039320969442e-17, -3.9630951242056644e-18, 1.9779107704573469e-15,
> > 7.8103367010164485e-15, -1.2351275675875924e-17, -8.6959160183013082e-15,
> > -1.6013333002787863e-18, 3.0110116065267407e-19, 3.7155867714997651e-16,
> > -0.12490087179942209, -2.2029631919898945e-16, -8.0869824211018481e-31,
> > -5.0586739111186695e-33, 2.0667225433033359e-19, 1.0141735013231455e-18,
> > 3.2845977347050714e-17, -2.9662734037779067e-19, 9.9595082331331027e-05,
> > -1.9982466114987023e-18, 2.6420113469249318e-18, 5.2415327077690342e-19,
> > 3.007533948014542e-18, -5.1128158961673558e-17, 6.7791246559883483e-19,
> > -3.042327665547861e-15, 1.9523738271941284e-19, -4.0556768902105205e-20,
> > -1.024372770865438e-17, -0.010638955366526863, 2.5786757042756338e-17,
> > 1.3546130005558746e-32, 1.7947249868721723e-33, 1.9228405363682382e-19,
> > -4.1209014001548407e-18, 1.5228291529275058e-17, 2.0443162438349059e-17,
> > -1.9572693345249383e-18, 0.0012709581028842473, -0.0018515998737074948,
> > 9.3128024867175073e-20, -5.1895788671618993e-17, 2.7373981615289174e-16,
> > 1.2812977711597223e-17, -2.2792319486263267e-15, 4.1599503721278813e-19,
> > -3.7733269771394201e-20, 4.16234419478765e-17, -1.5986158133468129e-16,
> > -0.016037670900735834, -5.0763064708173244e-33, -1.0176066166236902e-50,
> > -5.9296704797665404e-18, -8.0698801402861772e-17, 2.9646619457173492e-16,
> > -2.8879431119718227e-17, 2.9393253663483117e-18, -0.0018515998737074959,
> > 0.090335687736246548, 5.6849412731758135e-19, 8.8934458836007799e-17,
> > -4.1390858756690365e-16, 4.120323677410211e-16, 2.8000915545933503e-15,
> > 2.8094462743052983e-17, 1.1636214841356605e-18, 8.1510367497000071e-16,
> > -3.004558574117108e-15, 0.0061666744014873013, -2.5381532354086619e-33,
> > -2.7110175619881585e-34, -3.9720857469692968e-18,
> -2.2722246209861298e-17,
> > 6.5087631702810744e-18, -5.8833898464959171e-18, 6.0598974659958357e-19,
> > -6.4324667436451943e-19, -4.9865458549184915e-19, 0.010690736164778532,
> > -2.5717220776886191e-17, -2.9932234433250055e-17,
> -1.7586566364625091e-32,
> > -2.0572361612722435e-15, -4.1554892682395136e-18, 7.7947068522170098e-19,
> > 2.2950757715435305e-16, -6.8563401956907788e-17, 8.3986032618135276e-18,
> > -3.0929447793865389e-45, -9.1799940583213546e-47, 1.9205112318761867e-18,
> > -6.309535361758424e-17, -8.5794270550008683e-16, 1.9513428795022368e-15,
> > 6.975277535512248e-18, -6.2501857399599432e-17, 1.0008535047159303e-16,
> > -2.6690493301637561e-17, 0.11645557445978426, -5.0477656527863884e-17,
> > -1.9142598204346162e-30, -4.3865060631657549e-14,
> -9.2939338818469881e-18,
> > -3.7687560169835257e-19, 6.3729886582182645e-16, -1.2613712396302566e-14,
> > 7.8691203779893581e-16, -7.96410737362738e-44, -5.4321486152322576e-46,
> > -1.7940674523534164e-17, -8.157501463813461e-17, -3.3732156955561982e-15,
> > 7.677772572020035e-15, -4.7764170203097078e-17, 2.7500679201468897e-16,
> > -3.9378464227822661e-16, -1.9805344603340844e-17,
> -2.0590552280380157e-16,
> > 1.7227826719189774, -2.5152187296963038e-30, 1.7551764763918126e-14,
> > -5.4770085336579068e-18, 3.5206263799487875e-18, 8.2395392572605337e-16,
> > -4.1620820803543103e-14, -3.4715380153547894e-15,
> -3.0529078265571622e-43,
> > -2.8046838240348109e-45, -2.5479250226286821e-32, 4.121297503793449e-18,
> > 1.5823980389961397e-17, -1.235127567587614e-17, 1.0285610559852655e-18,
> > 1.0765062746617977e-17, 4.1501588000780158e-16, 5.2738511582000153e-33,
> > -1.6831260511527137e-30, -3.3541828527847257e-30, 3.7154414411814374,
> > 3.4040266595702152e-15, 5.1153310393987592e-18, 4.9999747986237238e-33,
> > -4.1627442819337855e-17, -1.3972969142201618e-16,
> -2.2035880670651683e-16,
> > -7.2138653746624827e-47, 2.5041819857453665e-47, -1.0640003589174214e-15,
> > -7.9766775252940816e-15, -1.86380928624219e-15, -2.0409565483972152e-14,
> > -4.5010528758464713e-16, -2.5476765314056885e-15, 5.7438239741680622e-15,
> > -5.1005913365528254e-15, -4.0052224901292309e-14, 2.0275682895130953e-14,
> > 8.001473599959262e-15, 4342.0489349328445, -2.2045103993340862e-15,
> > 2.087963708926218e-16, 8.0568968211307842e-14, 2.8167998366347968e-13,
> > 3.1749229442581484e-14, 8.9567352491809379e-43, -2.3670298646799273e-44,
> > 4.2781010071619158e-32, 2.5122662667827335e-17, -3.4359127754145909e-17,
> > -1.4706600926070093e-18, 1.876439844639638e-19, -9.1730442893545148e-20,
> > 2.8080363645722388e-17, -4.4705925750481468e-32, -7.809262716208445e-18,
> > -5.8577085238749991e-19, 5.1153310393987515e-18, -6.9351338840590931e-16,
> > 0.012093826986889086, -8.3952223993262708e-33, -2.5375314514710393e-16,
> > 3.424781603493142e-16, -4.3266783076648891e-18, 7.3570155156909413e-45,
> > 1.6171838713232998e-46, -0.00092251429799113343, -7.8240087156410085e-18,
> > 2.1083287560348406e-17, -8.9775735933099227e-18, 3.2531284806658474e-20,
> > -2.3137400470270888e-19, 7.166222774364368e-19, -1.6962655186338848e-19,
> > 2.5417429127255315e-18, -2.1888401697215364e-19, 7.492787666039381e-33,
> > 2.1373531604105888e-16, 5.2592866998595688e-19, 0.0053508323628379687,
> > 6.939324860007723e-17, -1.1345120732507663e-16, -1.3255448103671732e-18,
> > 1.0169078705557416e-16, -3.701528155571061e-19, -1.3050383760465465e-16,
> > -0.129172198695584, -1.5780816956514754e-14, -1.1071526119482569e-15,
> > -2.7949359117065728e-16, 9.8321930809926744e-17, 7.5470671367369062e-16,
> > 1.9622092961540385e-16, 7.9567529294094524e-16, 4.1183571472477233e-16,
> > -4.1627442819339772e-17, 1.2120022499676643e-13, -2.3248889366737827e-16,
> > 1.4240632385802796e-17, 1.3217752807216008, 2.0311176273844298e-13,
> > -1.9570716831733758e-15, 4.7962597568986823e-16, 2.9165264143720962e-17,
> > 2.5255547350017705e-16, -7.2487511301945062e-15, -0.44701112744185578,
> > -0.12490087179941219, -0.010638955366526433, -7.9593195538363592e-17,
> > -3.0600743713008676e-15, -1.8184344227119647e-16,
> -1.2723803109870239e-14,
> > -4.267580759582334e-14, -1.0227267431703487e-16, 3.9402946853735757e-13,
> > 3.1994002076770255e-16, -4.2061670894524917e-17, 7.4223786548689626e-14,
> > 7.0315216015826767, 2.2390669558248378e-15, -1.2244512793012972e-15,
> > -3.5253412228046429e-17, -1.2061427540284553e-18, 9.9899822637462946e-17,
> > -2.4550261028932595e-16, -2.5317628873304297e-16, 2.5199623506843722e-17,
> > -0.016037670900735834, 0.0061666744014872909, -1.3133261974683634e-18,
> > 6.5223470430906591e-16, -3.4512393153190672e-15, -2.4620056039335759e-16,
> > 2.8862767590185475e-14, -1.0881366439433239e-17, -3.9120910874012715e-18,
> > -1.2372606267786687e-15, 3.2639925872240233e-15, 0.30212469777249018,
> > 4.3378386520217752e-16, 1.1399322374260218e-17, 2.5273795605894831e-18,
> > 7.0372281397307745e-16, -1.4493840318281129e-15, 2.3909766400689098e-16,
> > -2.3621724602135697e-17, 4.1164300362570093e-18, -5.9299103903990776e-18,
> > -1.6549632537822852e-30, 2.0576427541669881e-29, 1.1492923590019121e-28,
> > -3.4572603795675457e-31, -2.0561228611852425e-29, 2.1007427249504219e-30,
> > 9.6204477501104748e-17, -6.7166997975872844e-15, 1.3428325376228686e-14,
> > 3.3362903423571782e-16, 3.2947112676731014, 6.6046445483993014e-18,
> > -4.0887487160499735e-20, -3.0968390611719402e-18,
> -6.0930441874738707e-19,
> > -3.1404422740276848e-18, -5.6584069489838701e-20, 1.6491214275041756e-19,
> > 6.3584544371811611e-22, 4.4921015596154276e-33, -3.1054789027110318e-31,
> > -1.1913370661394943e-30, 3.0677051354054613e-33, 5.8136140932799355e-30,
> > -4.0289214710854143e-33, 6.392220812489477e-19, 6.2007829741095621e-17,
> > 1.0190707777780152e-17, 9.6388161040420896e-18, -4.3748284667364627e-18,
> > 0.0054985968634936964), dim = c(19L, 19L), class = c("matrix"))
> >
> > If I run the following several times, given that the same seed is used,
> > I expect the numbers generated are the same:
> >
> > set.seed(1234)
> > MASS::mvrnorm(n = 5, mu = mu0, Sigma = sigma0)[1, 1:5]
> >
> > However, this is not the case:
> >
> >> set.seed(1234)
> >> MASS::mvrnorm(n = 5, mu = mu0, Sigma = sigma0)[1, 1:5]
> > [1] 0.4851605 0.5704446 1.6873036 2.7645014 -0.2020908
> >> set.seed(1234)
> >> MASS::mvrnorm(n = 5, mu = mu0, Sigma = sigma0)[1, 1:5]
> > [1] 0.4851605 0.5704446 1.6872818 2.7644787 -0.2023346
> >> set.seed(1234)
> >> MASS::mvrnorm(n = 5, mu = mu0, Sigma = sigma0)[1, 1:5]
> > [1] 0.4851605 0.5704446 1.6872818 2.7644787 -0.2023346
> >> set.seed(1234)
> >> MASS::mvrnorm(n = 5, mu = mu0, Sigma = sigma0)[1, 1:5]
> > [1] 0.4851605 0.5704446 1.6873036 2.7645014 -0.2020908
> >
> > It seems that the numbers in some columns alternate between
> > two versions.
> >
> > I understand that MKL can yield random numbers different from
> > those generated in the "standard R", due to the output of
> > eigen() used in MASS::mvrnorm(). However, should the numbers
> > generated identical if the same seed is used?
> >
> > This is the sessionInfo() output. I used a fresh installation of Ubuntu,
> > installed MKL, and then compiled R-devel.
> >
> >> sessionInfo()
> > R Under development (unstable) (2023-08-15 r84957)
> > Platform: x86_64-pc-linux-gnu
> > Running under: Ubuntu 22.04.3 LTS
> >
> > Matrix products: default
> > BLAS/LAPACK: /opt/intel/oneapi/mkl/2023.2.0/lib/intel64/libmkl_rt.so.2;
> > LAPACK version 3.10.1
> >
> > locale:
> > [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
> > [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
> > [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
> > [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
> > [9] LC_ADDRESS=C LC_TELEPHONE=C
> > [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
> >
> > time zone: Asia/Shanghai
> > tzcode source: system (glibc)
> >
> > attached base packages:
> > [1] stats graphics grDevices utils datasets methods base
> >
> > loaded via a namespace (and not attached):
> > [1] MASS_7.3-60.1 compiler_4.4.0
> >
> > I cannot create the same MKL R-devel machine used in CRAN checks.
> > Therefore, I am not 100% certain whether the cause of the error
> > is the same. Nevertheless, this phenomenon, non-identical random
> > numbers generated even with the same seed, can explain the
> > error I encountered.
> >
> > Regards,
> > Shu Fai
> >
> > ______________________________________________
> > R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
>
> ______________________________________________
> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
[[alternative HTML version deleted]]
More information about the R-help
mailing list