6
$\begingroup$

I was making a code change that required me to change some code that looked like

Abs[list]

to something that looked like

Map[Abs, list]

During testing I discovered that these two methods gave slightly different results. Some code to illustrate:

l1 = RandomComplex[{-1 - I, 1 + I}, 1000000];
Abs[l1] - Map[Abs, l1] // Norm

The output I get from running the code above in Mathematica 12.0.0.0 is on the order of 10^-14. My understanding was Abs "threading" over a list was the same as mapping Abs on each element of that list so I was expecting an output of 0. Can someone correct my understanding or otherwise explain what is the difference between these two methods?

$\endgroup$
8
  • 1
    $\begingroup$ Those results are identical within numerical accuracy. See Abs[l1] - Map[Abs, l1] // Chop // Total which returns exactly $0$. $\endgroup$
    – MarcoB
    Commented May 23, 2022 at 15:29
  • $\begingroup$ Maybe relevant mathematica.stackexchange.com/questions/76896/… $\endgroup$ Commented May 23, 2022 at 15:43
  • 3
    $\begingroup$ @MarcoB I do understand they are within numerical accuracy but I'm not sure why numerical accuracy should be a factor at all. It seemed to me that Abs[list] should perform the exact same operations as Abs[#]&/@list. The same operations on the same inputs should get you the same results, floating point error and all. $\endgroup$ Commented May 23, 2022 at 20:59
  • 1
    $\begingroup$ I repeatedly get 0. on a Mac M1 Max. Are you using an Intel CPU? Vectorized calls call the vectorized function in the MKL, which runs a different library code than a call on an individual number. There's a SystemOptions[] that controls when vectorized calls are made. Raise it to over 1000000 to check if it fixes the problem. (I ran into an occasional 1-bit difference some years ago in some function, but I can't recall the details.) $\endgroup$
    – Michael E2
    Commented May 24, 2022 at 4:37
  • 1
    $\begingroup$ Found it: It was Sin[]. See chat.stackexchange.com/transcript/message/39144034#39144034 and chat.stackexchange.com/transcript/2234?m=39169902#39169902 $\endgroup$
    – Michael E2
    Commented May 24, 2022 at 4:40

1 Answer 1

4
$\begingroup$

I observed similar behavior in a number of vectorized functions around five years ago. Since I figured it out to the extent possible, without knowing the details of the Math Kernel Library and how it takes advantage of the architecture of Intel CPUs, I didn't have a question to ask on site. Instead, as I sometimes do in such cases, I posted this arcane behavior in Chat, in case it is of interest to others. As I expected, it's usually not, another reason for not creating a new Q&A. But now, this question turns up, and the old Chat posts seem to be useful. V13 in my Mac M1 Max does not replicate the Intel behavior, which I can no longer check. Thus the first point is that the behavior in the OP is (probably) system-dependent.

The upshot of what I discovered about Intel machines was that different MKL functions are used for vectorized and nonvectorized math functions. The vectorized one must certainly take advantage of the pipelining and other advanced features of Intel CPUs to maximize memory management and throughput. Consequently, one imagines that the order of operations might be different. If so, an error of <1 ulp is conceivable and perhaps acceptable in Intel's design. All this is speculation but consistent with the evidence below. In the example I saved in Chat, which deals with the same behavior of Sin[] as the OP observes in Abs[], I found an example where there is a difference. And the difference turns out to be one of rounding the exact value up versus rounding it down, an error in each case of <1 ulp.

I left comments under the question, which the OP has had a chance to see. No one has responded, so I'm posting an answer (1) because I think it is probably the correct explanation and (2) to draw attention to it so that others might be motivated to check their systems. Here are the Chat postings:

Chat 1 (Aug 1, 2017 2:29 AM):

Vectorized results are not always identical with the mapped function:

foo = Subdivide[0., 2., 1100];
Sin@foo - (Sin /@ foo) // Abs // Max
Sin@foo - (Function[x, Sin[x], Listable]@foo) // Abs // Max

(* 5.55112*10^-17
   5.55112*10^-17  *)

Chat 2 (Aug 2, 2017 1:39 AM):

The difference above in the values of Sin[] has to do with a difference in how the values are computed when vectorized and not vectorized. Here's a single instance:

t1 = N[t0 = 1880891293519271/2251799813685248];
s1 = SetPrecision[Sin@t1, Infinity];
s2 = SetPrecision[First@Sin@ConstantArray[t1,257], Infinity];
s1 - Sin[t0] // N[#, 16] &
s2 - Sin[t0] // N[#, 16] &

(*  5.567563500914014*10^-17
   -5.534666745337552*10^-17  *)

The fraction t0 is the exact representation of the binary fraction for the floating-point t1. The vectorized result s2 is slightly better than the singleton-computed s1. For most inputs to Sin[], the two values are bitwise identical.

The magic number 257 comes from:

SystemOptions["ParallelOptions"]

{"ParallelOptions" -> {...,
   "VectorArithmeticThresholds" -> {{256, 2147483647}, {256,
      2147483647}, {256, 2147483647}}}}

I don't know what the numbers represent, but if you change the second 256, you change the threshold needed for the length of ConstantArray in the example.

$\endgroup$
2
  • 2
    $\begingroup$ Thanks for posting this as an answer. I did take a look at my system options in 12.0 and what I see is {"ParallelOptions" -> ..., "VectorVendorLengthThresholds" -> {128, 32, 32}}}. If I change my example to use 31 complex numbers, the difference I get is 0. With 32 values I get the floating point error. I believe you are right about differing MKL calls. Thanks very much for your investigation and sorry I did not get a chance to reply to your comments until now. $\endgroup$ Commented May 25, 2022 at 17:56
  • $\begingroup$ @ergodic_tortoise You're welcome. No problem about the late response. I was just explaining why I was posting a speculative answer. Thanks for letting me know it identified the issue. $\endgroup$
    – Michael E2
    Commented May 30, 2022 at 0:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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