在 Cython 中处理分数的最佳方法是什么？

Question

我有一个 np 数组，其中包含分数形式的元素以实现机器精度。我想应用像高斯消元这样的线性代数过程。这是我到目前为止的 Cython 代码（请注意，它只显示了获取上三角形形式的步骤，但实际上并没有引用它）。

在 Python 中生成的数据：

size = 5
foo = np.array([[Fc(v).limit_denominator(100) for v in r]
                for r in np.random.randn(size, size)])
identity = np.array([[Fc(v) for v in r] for r in np.identity(len(foo))])
m_id = np.concatenate([foo, identity], axis=1)

赛通：

%%cython
import numpy as np
cimport numpy as np
from quicktions import Fraction as Fc
cimport cython

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
def invert_gaussian4(np.ndarray matrix):

    cdef int matrix_size = matrix.shape[1] // 2
    cdef int c_i
    cdef int r_i
    cdef int swap

    for c_i in range(matrix_size - 1):
        swap = np.argmax(np.abs(matrix[c_i:, c_i])) + c_i
        matrix[[swap, c_i]] = matrix[[c_i, swap]]
        row = matrix[c_i, :] / matrix[c_i, c_i]
        for r_i in range(c_i+1, matrix_size):
            del_row = row * matrix[r_i, c_i]
            matrix[r_i, :] = matrix[r_i, :] - del_row

与 Python 相比，Cython 函数的性能并没有提高那么多。我已经认识到循环和分数元素中的 np 函数调用是减慢代码的原因。有关如何更好地优化此代码的任何建议？

Answer 1

您可以提高速度，也可以保持极高的精确度。你需要做出决定。

您在 cmets 中说过此数据用于“研究”，但未指定该字段。在相当多的研究领域中，产生的数据并不准确。相反，它们是您通过测量现实世界现象获得的近似值。我们说每个值都有多个有效数字。这些有效数字通过计算传播，最后，您应该四舍五入任何不重要的数字。

虽然浮点数学确实涉及中间舍入，但这种舍入通常会保留足够多的有效数字，最终结果不会受到影响。例如，使用 64 位双精度 IEEE 754 浮点值（默认情况下 Python 会这样做），您的有效数字有 53 个以 2 为底的有效数字，这等于以 10 为底的approx. 15 有效数字。如果您的实际数据例如，只有五个有效数字，那么对于大多数合理的操作，您不应该关心这种中间舍入。

如果这一段确实描述了你在做什么，那么你应该用标准浮点数替换你的分数对象。仅此一项就可以极大地加快您的计算速度。

如果这不能准确描述您在做什么（例如，因为您的领域是某种纯数学或其他理论学科），那么您可能别无选择，只能忍受缓慢。您可能会查看SymPy，它旨在执行这些领域倾向于关注的那种符号操作。

如果您使用的是真实数据，但它的有效数字超过 15 个，那么您应该知道 64 位整数只能达到 ~10^18。这意味着您的分数对象很可能是使用任意精度整数实现的，这确实非常慢。在这种情况下，您希望使用支持 128 位整数和/或浮点数的（超级计算）平台，并且您可能不想在 Python 中编码（阅读：预编译的 Python 二进制文件可能会或可能会这样的平台不存在，并且根据其标准一致性，它们可能会也可能不会自己编译；无论性能充其量是有问题的）。

最后，您不应该编写自己的高斯消除例程。相反，请使用numpy.linalg.solve。这可能会更快更精确。

Answer 2

此答案不会解决问题，但会解释您遇到问题的原因，并可能为您提供一些关于您需要妥协的线索。本质上你有两个问题：

1。您正在使用 Python 对象的 numpy 数组

如果 numpy 数组是由单个整数或浮点类型组成的数组，Cython 最适合它们。您的数组由任意 Python 对象组成（因此在 C 中，数组的每个元素都是指向存储 Python 对象的单独内存位置的指针）。这意味着每个索引操作都需要一个 incref 和 decref 用于引用计数。

这也意味着对元素的快速 C 级操作不可用。当您执行添加时，它必须在对象上对 __add__ 进行字典查找，然后调用该 Python 函数。同样，调用np.abs 对数组的每个元素执行__abs__ 的字典查找，然后调用Python 函数。类似地，调用 argmax 涉及在每个元素上查找比较运算符（可能是 __lt__）的字典。显然，这些东西很耗时，尤其是与纯 C 数组相比，其中每个操作都可能接近一个处理器指令。

解决此问题的一种方法是创建两个具有固定大小整数类型（理想情况下为 64 位）的数组来表示分数的分子和分母。 （注意：numpy 数组的类型由dtype 给出）。您需要自己在 Cython 中实现各种算术运算，但它是一个已知的数据类型这一事实会给您带来显着的加速（或者您可以使用 numpy 结构化数组将它们存储在一个数组中）。

但是：

2。 Fraction 类使用任意长度的整数

为了避免四舍五入，分数类使用 Python 的“长”整数，它可以容纳任意大的数字（达到内存的限制），但需要未知数量的内存这样做。

这样做的结果是不可能将它们分配到具有固定大小类型的连续数组中（这有利于索引速度）。第二个后果是不可能为它们生成“快速”的 C 代码 - 在每个操作期间，它需要计算当前的数字有多大，注意溢出，可能分配更多的内存。

最好的解决方案是接受您需要在某个地方舍入并使用固定大小的数据类型。第二种解决方案可能是有一个明确的“快速路径”和“慢速路径”。存储 4 个数组：分子小、分子大、分母小、分母大，小的为dtype=np.int64，大的为通用dtype=object。尽可能使用小的，并且只使用发生溢出的大的。如果您的大部分数据都是“小”的，这可能会更快。任何一种选择都需要自己编写大量算术代码。