高效的多精度数值数组答案

【问题标题】：Efficient multiple precision numerical arrays高效的多精度数值数组
【发布时间】：2015-08-01 17:48:52
【问题描述】：

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

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

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

要求“与本机浮点数一样高效”是没有意义的，但您可以要求它接近等效 C 代码（或者，如果不这样做，Java/C# 代码）的效率。特别是，高效的多精度数字数组意味着您可以进行矢量化操作，而不必在 Global Interpreter 中查找 __add__ 一百万次。

编辑：致亲密的选民：我的问题是关于将它们组合在一起的有效方式。可能重复中的答案特别指出幼稚的方法效率不高。

拥有一个 dtype=object 的 numpy 数组可能会产生误导，因为强大的 numpy 机制可以使标准 dtype 的操作超快，现在由默认对象的 python 运算符处理，这意味着速度将已经不在了

【问题讨论】：

可能相关：stackoverflow.com/questions/15307589/…
@cel 相关，但不相关。我的问题是关于效率，而答案提到这样做的效率低下那种方式。
和这个关系更密切：stackoverflow.com/q/26600471/1461210
@ali_m 这是另一个不做什么的例子，而这个问题问的是做什么。
是的，当然。我并不是说您的问题是重复的 - 我只是认为为了其他读者的利益，链接到 SO 上的其他相关问题很有用。

标签： python numpy mpmath gmpy

【解决方案1】：

免责声明：我维护gmpy2。以下测试是使用开发版本进行的。

a 和 b 是 1000 个元素列表，包含 250 位精度的伪随机 gmpy2.mpfr 值。测试执行两个列表的元素相乘。

第一个测试使用列表推导：

%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 开销，但随着精度的提高，有效节省的费用会减少。

【讨论】：

【解决方案2】：

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

浮点 numpy ndarray 中的每个元素占用相同数量的字节，因此该数组可以用第一个元素的内存地址、维度和连续之间的常规字节偏移量（或步幅）来表示数组元素。

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

相比之下，包含多个精度值的数组必须包含大小不等的元素，因为精度较高的值会比低精度值占用更多的字节。因此，多精度数组不能有规律地跨步，并且失去了上面提到的好处。

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

我怀疑这些性能问题可能是目前还没有提供您正在寻找的功能的 Python 库的主要原因。

【讨论】：

要求“与本机浮点数一样高效”是没有意义的，但您可以要求它接近等效 C 代码的效率（或者，如果不这样做，Java/C#代码）。特别是，高效的多精度数字数组意味着您可以进行矢量化操作，而不必在 Global Interpreter 中查找 __add__ 一百万次。
当然，我当然不是说不可能比当前情况下做得更好，即 GMP 标量单独包装在 Python 中。在您的问题中，您特别提到将 numpy 和 gmpy 结合起来作为提高性能的一种方式，所以我的回答主要是为了解释为什么我认为这不太可能奏效。
我的问题是关于如何比天真地结合 Numpy 和 gmpy 做得更好（即让 Numpy 将它们视为普通的 Python 对象）。另一种说法：将它们结合起来最有效的方法是什么？
不幸的是，没有比使用 np.object 数组（本质上将 gmpy 标量视为单独的 Python 对象）更有效的方式来组合 numpy 和 gmpy。为了更好地使用 numpy，有必要设计一种全新类型的数组容器（想想scipy.sparse 矩阵与np.ndarray）。你还没有真正描述你的用例是什么，但也许你最好的选择可能是用 C/C++/Cython 编写自己的扩展，直接调用 GMP 库函数。
我维护gmpy2。我已经研究过添加一些gmpy2 函数的矢量化版本。我使用gmpy2.sin、map(sin, list) 和gmpy2.sin(list) 比较了列表理解，并没有太大差异（大约 10%）。我没有尝试过使用nd.array 或基本的乘法/加法函数。如果你想看测试代码，请告诉我。

【解决方案3】：

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

（如果您想查看 dtype 的外观，您可能已经取消了使用 Numpy 支持编译它的相关行的注释；它应该可以用于一目了然，但尚未实现任何功能；下一个应该在 9 月或 10 月发布。）

【讨论】：