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

Question

Numpy非常适合通过矢量化来加速代码。然而，当执行例如隐式迭代方程求解方法(如高斯-塞德尔方法(时，则阵列的元素n+1取决于n的第一元素，从而阻止矢量化。

我已经实现了这样的高斯塞德尔：

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循环，但这是一种加快速度的方法吗？

Answer 1

将if测试转移到以下循环之外可能会有所帮助

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]

如果可以避免测试的话。

另一个需要尝试的东西是numba包，它可以在没有外部编译步骤开销的情况下动态编译代码。