并行计算核密度估计时结果与串行不一致的问题咨询

我希望在网格上并行执行核密度估计（KDE），但遇到缩放问题，多线程计算结果无法和串行计算保持一致。以下用Scipy官方文档的示例复现该问题：

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

def measure(n):
    "Measurement model, return two coupled measurements."
    m1 = np.random.normal(size=n)
    m2 = np.random.normal(scale=0.5, size=n)
    return m1+m2, m1-m2

m1, m2 = measure(2000)
xmin = m1.min()
xmax = m1.max()
ymin = m2.min()
ymax = m2.max()

kde_bandwith = 0.5

X, Y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]
positions = np.vstack([X.ravel(), Y.ravel()])
values = np.vstack([m1, m2])
kernel = stats.gaussian_kde(values, bw_method=kde_bandwith)
Z = np.reshape(kernel.evaluate(positions).T, X.shape)

我将最后两行的计算拆分到n_split个子数据集上执行（实际会在不同线程运行，此处用循环简化）：

## split
n_split = 10
values_split = np.array_split(values, n_split, axis=1)

results = []
for v in values_split:
    kernel_n = stats.gaussian_kde(v, bw_method=kde_bandwith)  
    Z_n = np.reshape(kernel_n.evaluate(positions).T, X.shape)
    results.append(Z_n)

Z2 = np.sum(results, axis=0)

通过将每个子集的Z_n求和得到合并结果Z2，绘图对比结果：

• 左图：串行计算的原始结果
• 中图：子集合并后的结果
• 右图：两者的差值

注意：KDE的幅度看起来是n_split倍，但即使简单用1/n_split或len(values_split)/len(values)（本例中两者相等）进行归一化，结果仍不完全一致，相关对比图如下：

内容的提问来源于stack exchange，提问作者Philippe Miron