Python验证序列循环排列引理时浮点数比较异常求助

背景 #

我正在验证一个数学引理：给定任意实数序列$a_1,…,a_n$，记其累积和为$s_1,…,s_n$，则存在唯一的循环排列$\sigma$，使得对所有$k$，$s_k(\sigma) ≤ \frac{s_n}{n} \cdot k$（即部分和序列位于连接$(0,0)$与$(n,s_n)$的直线下方）。

我编写了Python代码遍历所有循环排列（排除初始排列）来寻找符合条件的排列并绘图，但遇到了异常情况。

异常现象 #

• 使用np.random.randn(9)生成随机浮点数数组时，程序有时会判定不存在符合条件的排列，触发else分支打印序列。
• 但把此时输出的序列手动输入代码，却能正常找到并绘制目标排列。
• 使用np.random.randint(-10,10,9)生成整数数组时，程序始终正常运行。

代码与出错序列 #

验证代码 #

import numpy as np
import matplotlib.pyplot as plt

randarr = np.random.randn(9)

tot_arr = np.array([0, *randarr])

init_cumsum = np.cumsum(tot_arr)
mean_array = np.array([ n * init_cumsum[-1]/(len(init_cumsum) - 1)\
                      for n in range(len(init_cumsum)) ])

def permuted(k):
    return np.cumsum([0, *np.roll(randarr, k)])


fig, ax = plt.subplots(figsize=(5, 3))
ax.plot(np.arange(0, 10, 1), init_cumsum)
ax.plot(np.arange(0, 10, 1), mean_array, '--')

for k in range(1, 10):
    if (mean_array >= permuted(k)).all() :
            ax.plot(np.arange(0, 10, 1), permuted(k), label="correct permutation")
            plt.legend()
            break
            

else:
    print(tot_arr)

plt.show()

出错示例序列 #

-0.12254795 -1.18301478  1.9911016  -0.67975716 -0.53877141 0.49170832 -1.08431192  1.03930656 -0.14305318
-0.98415686 -0.19698914  1.25619013 -0.31378433  1.45470612 0.04794845 -2.12723127 -0.69520889  0.28457194
2.11822629 -1.2222578   0.48015473 -0.57059476  0.76872146 0.44806995 -0.76466846  0.404127   -0.35847348
-0.04126626 -2.0829324  -0.5147483   0.13324286 -2.08391787 0.00566639  2.07952122  0.83424157  0.71345859

原因分析 #

问题根源是浮点数精度误差：

1. 浮点数在累积和、均值计算过程中会产生微小的舍入误差，导致循环排列后的部分和序列permuted(k)与mean_array的比较中，出现本应相等的位置，实际计算结果是permuted(k)略大于mean_array的情况。
2. 代码中使用(mean_array >= permuted(k)).all()进行严格的逐元素比较，只要有一个位置不满足，就判定该排列不符合条件，从而错过正确的排列。
3. 手动输入序列时，打印出来的数值是经过截断的近似值，重新计算时的累积和误差可能刚好不会触发这个不满足的情况；而整数序列的计算是精确的，不会有精度问题，所以始终正常运行。

修复方案 #

将严格的逐元素比较替换为允许微小误差的近似比较，使用np.allclose函数：

# 替换原有的if判断
if np.allclose(permuted(k), mean_array, atol=1e-8) or (permuted(k) <= mean_array + 1e-8).all():

或者直接给比较加一个微小的容差：

if (mean_array + 1e-8 >= permuted(k)).all():

这样就能忽略浮点数计算带来的微小误差，正确找到符合引理条件的循环排列。

内容的提问来源于stack exchange，提问作者Daigaku no Baku