## Multiplication non-deterministic?

3

2

I have two vectors of 134 elements each ($mu$, and $gt$). $mu$ contains Integers, and $gt$ contains machine precision Reals. I execute the following simple expression multiple times without changing either mu or gt: $$(mu/2*gt).gt$$ I will get one of two different results: $88474.52216839303$ or $88474.52216839301$ (they differ in the last digit). What's up with that? I have never seen Mathematica behave non-deterministically on simple mathematical expressions before (I'm using version 9.0.1). Could it be performing the dot product or the multiplication in parallel or something? I really need deterministic results (at least on the same computer). I'm using Windows 7 on a Dell laptop.

Update 2. More fun (but apparently only when n = 134).

n = 134;


Manually evaluate the following expression over and over again in a notebook. I get a different output for Hash[a.a] about every 5 or 6 evaluations (I realize that I'm using RandomReal, but I'm resetting the seed every time so $a$ always receives the same value).

SeedRandom[1]
a = RandomReal[{0, 1}, n];
Hash[a]
Hash[a.a]


Update 1. Here's the output of two supposedly identical computations:

InputForm[(mu/2 gt).gt]
88474.52216839301

InputForm@Sum[(mu[[i]]/2 gt[[i]]) gt[[i]], {i, Length@mu}]
88474.52216839303


Sometimes the first expression comes out identical to the second expression.

The two vectors are:

InputForm[mu];
{1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000,
1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000,
1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000,
1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000,
1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000,
1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000,
1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000,
1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000,
1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000,
1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000,
1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000,
1000, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10};

InputForm[gt]
{-0.004653655618380714, -0.009080939192239956, -0.031034220835509885,
0.06650237015337201, -0.009447749078904444, 0.03336916426793746,
-0.4631539984577655, -6.265611574451846, 0.5050454470606143,
-0.19087962150762205, -4.3828201756709735, 0.49876272299955815,
0.3861679596859371, -6.1635962937354645, -0.15452549964472695,
0.3473530604725559, -4.199256945736951, -0.040039918046105516,
0.004474673856913913, -0.08203350223995054, 0.00755080417433085,
0.19742429370731282, -0.058492983677093546, -0.07179532877025352,
-0.06654467292071248, 0.053278032391796154, 0.021824097049037344,
-0.15494663566691225, 0.0490691048534343, 0.06064592126063756,
-0.07745582064698114, -0.00746967588173042, 0.03559995124065775,
-0.013295758759910827, -0.0021771809232924225, 0.014751625280327196,
-0.012490576622492892, -0.004059719351428424, 0.02260878725186556,
0.00906537308368245, 0.00005640969962929232, 0.013322793809346312,
0.01939688322252025, 0.0012921413237499124, 0.005353410478874665,
0.013744515181619552, -0.0008260952031691714, 0.0019523689330869926,
0.0036484573285596087, 0.0002894048077433703, 0.0009607045972313589,
-0.0005226991236028443, -0.002322046639359741, -0.024682514037441098,
0.027220500404569047, 0.0545967335062756, 0.05642841641401902,
0.004367474285221214, -0.0005337849541497026, -0.0002226867162167549,
0.00019518810603202322, 4.710483244170917*^-6,
-0.00016670488036184576, -0.00011284285936618362,
-0.00008617428801798452, -0.00009997965667755927,
-0.00011813749998897688, -0.00011200191449290103,
-0.00003504307476243085, -0.00007170970464401294,
8.314045093271716*^-6, -0.000012108637423846602,
-0.00008318583029878757, 6.872184070694232*^-6,
0.000039598867001222615, -0.00030723457677017674,
0.000021126291712914025, 0.00017867733864129764, -0.00132238312878738,
0.00008352828516677846, 0.00016343481838493323,
-0.0009495187710002947, 0.000014501494743246468,
0.00010295794136475078, -0.0007041704456466069,
-0.000038074551413580515, 0.000024071562009758685,
-0.000368866036143653, -0.000035157347339220824,
-7.738001399044264*^-6, -0.00012744527523419168,
-0.00001767740166705445, -0.00006515075020033605,
-0.0005012315885292151, -0.00009246376382040004,
0.000035397133843528654, 0.00010514304758398063, 0.000164070573737746,
0.00033034132631488825, 0.0002563940375876836, 0.0001581254676494522,
0.00010340418143852661, 0.0001271214003914325, 0.000144619937452181,
-6.445829617805765*^-8, -7.918633438308738*^-8,
-5.331850464748234*^-10, -5.395049641407357*^-11,
-9.692675239369703*^-10, -2.3443641234029144*^-10,
-2.741707881994002*^-9, -1.8602797897764294*^-9,
1.852265547443208*^-10, -8.268466237349803*^-9,
-2.873207919632894*^-10, -8.965460213556384*^-9,
-2.7681125512775745*^-9, -6.542104162408522*^-9,
-4.0174693373664825*^-9, -3.665420456657473*^-9,
-2.843221227293432*^-9, -1.5467319689177169*^-9, 0.12984908695319053,
0.3495235744679681, 0.31365238179137667, 0.5566921439809294,
0.09068947055843457, 78.52618761194245, 1.462920938774925,
0.13819238950337862, 0.24562317270333922, 0.28816819236661184,
0.02950510971255533, 0.19955099011468305}


1Can you post values for your vectors? If it is too long for this post, then you can put it for instance on pastebin. It would be really helpful if we had the vectors to try this. – halirutan – 2016-06-17T20:11:25.420

I updated the original post with the vectors. Thanks for taking a look. – Eric Parker – 2016-06-17T20:49:03.597

1I can not reproduce the anomaly with v10 (get ..301 every time). How many re-evals do you need to see to expect the different result? – george2079 – 2016-06-17T20:58:58.980

It does sometimes get "stuck" on one value or another for a long period of time. But after a kernel relaunch I can get it to return different results again. I should mention that I have a somewhat large application that produces these two vectors. I assign to global variables from the local variables within the app at the point where this dot product expression is calculated. Then I manually evaluate the expression using the global variables in a separate notebook. It could be that something in my app is affecting expression evaluation, but I'm not overloading any Mma functions. – Eric Parker – 2016-06-17T21:11:51.567

@george2079 You get ...301? Can you try evaluating the second expression in the post (the Sum[...])? – Eric Parker – 2016-06-17T21:15:07.317

i get the same results as in the question for the two different expressions. Maybe you can do an extended precision calculation InputForm[(mu/2 # ). # &@SetPrecision[gt, 20]] -> ...3019158 – george2079 – 2016-06-17T21:43:49.913

2Is the fluctuation still present if you do SetSystemOptions["MKLThreads" -> 1]; ? – ilian – 2016-06-17T22:04:53.043

@george2079 I further simplified my example to just a simple Dot product of a random vector. Just click and evaluate the expression over and over again. See if you get different output in ver 10 (see Update 2 in the post) – Eric Parker – 2016-06-18T00:55:31.450

6

It's not a bug and it's not so uncommon. For an explanation have a look here. This and some related issues also appear in this MathGroup thread.

Also relevant: 1 2.

Thanks for the references. I understand about machine numerics (e.g. 0.1). My issue is with determinism. I believe that it is a bug that identical input can produce different results on the same platform, compiler, version, etc. I think that I am probably running up against the MKL ddot alignment issue that you documented. Here's a simple fix: just always allocate numeric arrays aligned to 16-byte boundaries inside Mma (most compilers have alignment specifiers). Anyway, thanks for explaining what is happening. – Eric Parker – 2016-06-18T13:53:29.250

3

They are not identical computations. With the first form,

(mu/2 gt).gt


Mathematica can take advantage of vector arithmetic, usually going through specialized routines like LAPACK. The second form,

Sum[(mu[[i]]/2 gt[[i]]) gt[[i]], {i, Length@mu}]


however, will usually be calculated term by term because there is a possibility that the input can change from term to term. The fact that they differ only in the last decimal place is surprising as various numerical errors can creep into summations if they are not carefully accounted for.

But even if it is true that the vector version might use a different code path, wouldn't you expect it to always use the same code path for the same inputs? I'm seeing different results randomly. – Eric Parker – 2016-06-17T20:47:07.347

2You're seeing variation in the last digit. That isn't non-determinism, that's noise in the floating point arithmetic. Ignore it. – rcollyer – 2016-06-17T20:59:45.453

3Floating point addition is only approximately associative, so your result will depend on the order of the summation. I suspect Dot recursively divides the array in halves, summing the halves (more accurate) while Sum does it term by term. – John Doty – 2016-06-17T21:36:13.497

1I am slightly amused that one is agonizing over an ulp. Nevertheless, one might also want to look at the result of using Total[] with "CompensatedSummation". – J. M.'s ennui – 2016-06-17T23:17:44.457