numpy与MATLAB归一化求和差异原因及Python适配MATLAB方法

问题现象 #

在Python 3.8 + NumPy中对形如$p_k := 2^{-k2}$（$k \in {1,\dots,n}$，$n=8$）的概率分布做归一化，使用等价于MATLAB num2hex的实现方法后，发现：

• 未归一化的分布在Python和MATLAB R2020a中完全一致
• 归一化后的单个分布值十六进制表示完全相同
• 归一化后分布的求和结果存在差异：Python输出0x3ff0000000000000（精确1.0），MATLAB输出0x3fefffffffffffff（比1.0小一个ULP的双精度值）
• 当$n<8$时，两者求和结果完全一致

Python代码及结果 #

import numpy as np
import struct # 实现MATLAB num2hex等价功能

n = 8
unnormalizedPDF = np.array([2**-(k**2) for k in range(1,n+1)])
# MATLAB num2hex等价实现
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

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'

差异原因 #

核心差异来自于浮点数求和的精度处理策略不同：

• MATLAB的sum函数采用逐元素递推累加，保留累加过程中产生的所有舍入误差，不会对最终结果做额外修正。当n=8时，归一化后的极小值累加带来的微小误差会被保留，最终结果略小于1.0。
• NumPy的np.sum内部使用了优化的累加算法（如利用更高精度的临时累加器，或对接近整数的结果做舍入修正），当求和结果非常接近1.0时，会直接将其舍入为精确的双精度1.0，从而消除了微小的误差。

n<7时，归一化后的所有值求和误差极小，不足以触发两种求和策略的差异，因此结果一致。

对齐MATLAB结果的解决方法 #

要让Python的求和结果与MATLAB一致，需要模拟MATLAB的逐元素递推累加逻辑，避免NumPy的求和优化。可以手动实现循环累加：

import numpy as np
import struct

n = 8
unnormalizedPDF = np.array([2**-(k**2) for k in range(1,n+1)])
num2hex = lambda x : hex(struct.unpack('!q', struct.pack('!d',x))[0])

# 先计算归一化后的数组
normalized_pdf = unnormalizedPDF / np.sum(unnormalizedPDF)
# 手动递推累加，模拟MATLAB sum行为
matlab_like_sum = 0.0
for val in normalized_pdf:
    matlab_like_sum += val

print(num2hex(matlab_like_sum))
# 输出：0x3fefffffffffffff，与MATLAB结果一致

该方法通过逐次将每个归一化值加到累加器中，完全复刻了MATLAB sum的递推逻辑，保留了累加过程中的舍入误差，从而得到和MATLAB完全一致的结果。

内容的提问来源于stack exchange，提问作者S Huntsman