[R] MASS::mvrnorm() on MKL may produce different numbers even when the seed is the same?

[Shu Fai Cheung]

Thanks a lot! I am not aware of this behavior! I will take this into
account in my work.

Regards,
Shu Fai

On Thu, Aug 17, 2023 at 10:31 PM Bill Dunlap <williamwdunlap using gmail.com> wrote:
>
> MKL's results can depend on the number of threads running and perhaps other things.  They blame it on the non-associativity of floating point arithmetic.  This article gives a way to make results repeatable:
>
> https://www.intel.com/content/www/us/en/developer/articles/technical/introduction-to-the-conditional-numerical-reproducibility-cnr.html
>
> -Bill
>
> On Wed, Aug 16, 2023 at 8:11 PM Shu Fai Cheung <shufai.cheung using gmail.com> 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.