高效的多精度数值数组

Efficient multiple precision numerical arrays

Numpy 是一个用于高效数值数组的库。

mpmath 在 gmpy 的支持下,是一个用于高效多精度数字的库。

如何有效地将它们组合在一起?还是仅使用带有 mpmath 数字的 Numpy 数组就已经很有效了?

要求 "as efficient as native floats" 没有意义,但您可以要求它接近等效 C 代码的效率(或者,如果做不到,Java/C# 代码) .特别是,高效的多精度数字数组意味着您可以执行矢量化操作,而不必在全局解释器中查找 __add__ 一百万次。

编辑:致最终投票者:我的问题是关于高效 将它们组合在一起的方法。 possible duplicate 中的答案特别指出了 naive approach 的效率不高。

Having a numpy array of dtype=object can be a liitle misleading, because the powerful numpy machinery that makes operations with the standard dtypes super fast, is now taken care of by the default object's python operators, which means that the speed will not be there anymore

据我所知,目前还没有 Python 库支持对多精度值进行矢量化数组操作。不幸的是,没有特别有效的方法在 numpy ndarray 中使用多个精度值,而且这种情况极不可能发生,因为多个精度值与 numpy 的基本数组模型不兼容。

浮点数 numpy ndarray 中的每个元素占用相同的字节数,因此数组可以根据第一个元素的内存地址、维度和常规字节偏移量(或跨度)来表示在连续的数组元素之间。

此方案具有显着的性能优势 - 相邻的数组元素位于相邻的内存地址,因此顺序 reads/writes 数组受益于更好的引用局部性。跨步对于可用性也非常重要,因为它允许您在不在内存中创建新副本的情况下对同一数组的视图进行操作。当您执行 x[::2] 时,您实际上只是将数组第一个轴上的步幅加倍,这样您就可以处理所有其他元素。

相比之下,包含多个精度值的数组必须包含大小不等的元素,因为较高精度值会比低精度值占用更多字节。因此,多精度数组不能定期跨越,并且失去了上述好处。

除了构造数组的问题之外,即使是多精度标量上的简单算术也可能比浮点标量慢得多。几乎所有的现代处理器都有专门的浮点单元,而多精度运算必须在软件而不是硬件中实现。

我怀疑这些性能问题可能是 Python 库无法提供您正在寻找的功能的主要原因。

免责声明:我维护 gmpy2。使用开发版本执行了以下测试。

ab 是包含 250 位精度的伪随机 gmpy2.mpfr 值的 1000 个元素列表。该测试执行两个列表的逐元素乘法。

第一个测试使用列表理解:

%timeit [x*y for x,y in zip(a,b)]
1000 loops, best of 3: 322 µs per loop

第二个测试使用map函数执行循环:

%timeit list(map(gmpy2.mul, a, b))
1000 loops, best of 3: 299 µs per loop

第三个测试是列表理解的 C 实现:

%timeit vector2(a,b)
1000 loops, best of 3: 243 µs per loop

在第三次尝试中,vector2 试图成为一个行为良好的 Python 函数。使用 gmpy2 的类型转换规则处理数字类型,完成错误检查等。检查上下文设置,根据要求创建次正规数,如果需要则引发异常等。如果忽略所有Python 增强并假设所有值都已经 gmpy2.mpfr,我能够在第四次尝试时缩短时间:

%timeit vector2(a,b)
10000 loops, best of 3: 200 µs per loop

第四个版本没有进行足够的错误检查以供一般使用,但第三次和第四次尝试之间的版本可能是可行的。

可以减少 Python 开销,但随着精度的提高,有效节省会减少。

当前的项目是 qd,它将能够利用其值在内存中的固定大小,将高精度数字嵌入到 Numpy 数组中。现在,该类型可用于 Numpy,但还不能作为 dtype;但是,您已经可以将它与对象 dtype 一起使用。

(如果您想查看 dtype 的外观,您可能已经取消注释相关行以使用 Numpy 支持对其进行编译;它应该可以一目了然,但没有任何功能尚未实施;下一个版本应该在 9 月或 10 月。)