1

I tried to normalize a probability distribution of the form $p_k := 2^{-k^2}$ for $k \in {1,\dots,n}$ for $n = 8$ in numpy/Python 3.8 along the following lines, using an equivalent of MATLAB's num2hex a la C++ / Python Equivalent of Matlab's num2hex. The sums of the normalized distributions differ in Python and MATLAB R2020a.

If $n < 8$, there is no discrepancy.

What is going on, and how can I force Python to produce the same result as MATLAB for $n > 7$? It's hard for me to tell which of these is IEEE 754 compliant (maybe both, with a discrepancy in grouping that affects a carry[?]) but I need the MATLAB behavior.

I note that there are still discrepancies in rounding between numpy and MATLAB per Differences between Matlab and Numpy and Python's `round` function (which I verified myself) but not sure this has any bearing.

import numpy as np
import struct # for MATLAB num2hex equivalent below
n = 8
unnormalizedPDF = np.array([2**-(k**2) for k in range(1,n+1)])
# MATLAB num2hex equivalent a la https://stackoverflow.com/questions/24790722/
num2hex = lambda x : hex(struct.unpack('!q', struct.pack('!d',x))[0])
hexPDF = [num2hex(unnormalizedPDF[k]/np.sum(unnormalizedPDF)) for k in range(0,n)]
print(hexPDF)
# ['0x3fec5862805436a4', 
#  '0x3fbc5862805436a4', 
#  '0x3f6c5862805436a4', 
#  '0x3efc5862805436a4', 
#  '0x3e6c5862805436a4', 
#  '0x3dbc5862805436a4', 
#  '0x3cec5862805436a4', 
#  '0x3bfc5862805436a4']
hexPDFSum = num2hex(np.sum(unnormalizedPDF/np.sum(unnormalizedPDF)))
print(hexPDFSum)
# 0x3ff0000000000000

Here is the equivalent in MATLAB:

n = 8;
unnormalizedPDF = 2.^-((1:n).^2);
num2hex(unnormalizedPDF/sum(unnormalizedPDF))
% ans =
% 
%   8×16 char array
% 
%     '3fec5862805436a4'
%     '3fbc5862805436a4'
%     '3f6c5862805436a4'
%     '3efc5862805436a4'
%     '3e6c5862805436a4'
%     '3dbc5862805436a4'
%     '3cec5862805436a4'
%     '3bfc5862805436a4'
num2hex(sum(unnormalizedPDF/sum(unnormalizedPDF)))
% ans =
% 
%     '3fefffffffffffff'

Note that the unnormalized distributions are exactly the same, but the sums of their normalizations differ by a single bit. If I use $n = 7$, everything agrees (both give 0x3fefffffffffffff), and both give the same results for $n < 7$ as well.

S Huntsman
  • 123
  • 5
  • 1
    If your code depends on specific floating-point rounding errors, you’re doing it wrong. This is not only about differences in computation order, but also about differences in floating-point operations in different chips. Don’t try to get bit-level equality between programs, make your algorithm robust against rounding errors. – Cris Luengo Dec 02 '22 at 18:46
  • @CrisLuengo -- I am doing it right, but I agree with you in general. The MATLAB version (with direct summation) is correct. I am computing bounds on an (at times astronomically large) expected value and have proved several lemmas underlying these bounds. I wasn't looking for bit-level equality between programs but the Python np.sum() stuff was giving me a denominator of 0. If Python behaved properly and called the result Inf I probably would never have noticed any of this in the first place, but instead I got NaNs for my bounds that were totally unacceptable. – S Huntsman Dec 02 '22 at 19:09
  • PS- Direct summation is numerically ideal because the thing I'm summing is already sorted in descending order to make the numerics nice. – S Huntsman Dec 02 '22 at 19:10

1 Answers1

3

According to the manual, numpy.sum uses pairwise summation to get more precision. Another common algorithm is Kahan summation.

Anyway, I wouldn't count too much on Numpy and MATLAB giving the same result up to the last bit, as there might me operation reordering if computations are made in parallel. See this for the kind of problem that can arise.

However, we can cheat a little bit and force Python to do the sum without the extra precision:

hexPDFSum = num2hex(np.sum(np.hstack((np.reshape(unnormalizedPDF / np.sum(unnormalizedPDF), (n, 1)), np.zeros((n, 1)))), 0)[0])
hexPDFSum
'0x3fefffffffffffff'
  • Thanks for this! I can't tell what MATLAB uses, though it seems likely that it's direct summation in this case as suggested by https://www.mathworks.com/matlabcentral/answers/550-compensated-summation-in-sum. Summation is apparently done with more finesse for large arrays as https://blogs.mathworks.com/loren/2021/09/07/a-faster-and-more-accurate-sum/ spells out. – S Huntsman Dec 02 '22 at 15:54
  • I got Python to agree with MATLAB by using sum instead of np.sum as well. Thanks again! – S Huntsman Dec 02 '22 at 15:58
  • @SHuntsman I added a way to get the same result: by forcing the sum on a dimension that is not the leading one. Probably not efficient, but if you need it, you have it. –  Dec 02 '22 at 15:58