如何使用 numpy 为隐式方法快速编写循环答案

【问题标题】：How to write fast for loops with numpy for implicit methods如何使用 numpy 为隐式方法快速编写循环
【发布时间】：2020-02-25 12:02:05
【问题描述】：

Numpy 非常适合通过矢量化加速代码。但是，在执行例如 Gauss-Seidel method 之类的隐式迭代方程求解方法时，数组的元素 n+1 依赖于 n 的第一个元素，从而阻止了矢量化。

我已经这样实现了 gauss-seidel：

def gaussSeidel(A, b, x0, tol = 1e-10):
    n = len(A)
    x = np.copy(x0)
    k = 0

        while (True):
        print(k)
        for i in range(n):
            x[i] = (b[i] - np.delete(A[i], i) @ np.delete(x, i)) / A[i, i]

        if max(abs(x - x0)) < tol:
            return x

        x0 = np.copy(x)
        k += 1

这种方式比显式版本慢，例如 Jacobi 方法，

def jacobi(A, b, x0, tol = 1e-10):
    n = len(A)
    M = np.copy(A)
    DInv = np.zeros((n, n))
    k = 0
    for i in range(n):
        DInv[i, i] = 1 / M[i, i]
        M[i, i] = 0

    while (True):
        x = DInv @ (b - M @ x0)
        if max(abs(x - x0))  < tol:
            return x
        x0 = x
        k += 1
        print(max(abs(x - x0)))

据我所知，这里不可能消除 for 循环，但这是一种加快速度的方法吗？

【问题讨论】：

您可以对内部 for 循环进行矢量化...这实际上是一个点积，省略了 i-th 元素...所以使用全点积并对其进行调整相应地
谢谢。我已经编辑了代码，它对一些人有所帮助，但主要问题仍然存在
能否添加与您相关的示例数据。 A,b,x0的数组形状是什么？
numpy 向量化的本质是将循环移动到已编译的代码 - 在整个数组方法中。
之前我建议使用点积并调整以删除内部循环...使用numpy.delete 这样做是不是的好方法，因为它会导致一份副本。

标签： python numpy numeric

【解决方案1】：

可能会有所帮助的事情是将 if-test 移到后续循环之外

        s = b[i]
        for j in range(n):
            if i == j:
                continue
            s -= A[i, j] * x[j]

如

        s = b[i] + A[i, i] * x[i]
        for j in range(n):
            s -= A[i, j] * x[j]

这样可以避免 n 个 if 测试。

要尝试的另一件事是numba 包，它可以即时编译您的代码，而无需外部编译步骤的开销。

【讨论】：