python - 为什么 Gauss-Jacobi 方法的特定 numpy 实现会显着减少迭代次数？

Question

在 python 中实现 Gauss Jacobi 算法时，我发现两种不同的实现需要截然不同的迭代次数才能收敛。

第一个实现是我最初想出的

import numpy as np
def GaussJacobi(A, b, x, x_solution, tol):
    k = 0
    N = A.shape[0]
    D = np.diag(A)
    R = A-np.diagflat(D);
    while(checkTol(tol, x, x_solution)):
        x_new = np.zeros(N, dtype=np.double) #x(k+1)
        for i in range(N):
            aii = D[i]
            bi = b[i]
            s = np.dot(R[i], x)
            x_n[i] = (1/aii)*(bi - s)
        x = x_new
        k+=1
        print('x(%d) =' % k, x)
    return k

第二种实现基于this article.

def GaussJacobi(A, b, x, x_solution, tol):
    k = 0
    N = A.shape[0]
    D = np.diag(A)
    R = A-np.diagflat(D);
    while(checkTol(tol, x, x_solution)):
        for i in range(N):
            x = (b - np.dot(R, x)) / D
        k+=1
        print('x(%d) =' % k, x)
    return k

解决下列问题时

A = [ 4, -1,  0, -1,  0,  0]
    [-1,  4, -1,  0, -1,  0]
    [ 0, -1,  4,  0,  0, -1]
    [-1,  0,  0,  4, -1,  0]
    [0,  -1,  0, -1,  4, -1]
    [0,   0, -1,  0, -1,  4] 

b = [2, 1, 2, 2, 1, 2]

x_solution =[1, 1, 1, 1, 1, 1]

x0 = [0, 0, 0, 0, 0, 0]

第一个实现需要 37 次迭代才能收敛，误差为 1e-8，而第二个实现只需要 7 次迭代即可收敛。

是什么让第二个实现比第一个实现快得多？

编辑:

我已经实现了另外两种方法，Gauss-Seidel 方法和 SOR 方法。这两种方法的实现方式与我最初的慢速 Gauss-Jacobi 方法类似。

我对 100 个 NxN 对角占优矩阵进行了随机测试，每个 N = 4...20 以获得收敛前的平均迭代次数。

  N    Gauss-Jacobi    Gauss-Jacobi Fast    Gauss Seidel    SOR -- w=1.5
---  --------------  -------------------  --------------  --------------
  4           40.96                17.04         40.6804         40.9204
  5           49.11                17.25         48.7489         48.9389
  6           56.11                16.04         55.6789         55.9089
  7           70.26                18            69.6774         70.0074
  8           76.4                 16.54         75.756          76.236
  9           83.56                17.03         82.8344         83.1044
 10           92.33                16.24         91.5267         91.7267
 11           98.02                16.59         97.1598         97.4598
 12          107.39                15.98        106.436         106.756
 13          123.48                17.75        122.375         122.655
 14          125.07                16.04        123.949         124.239
 15          132.41                16.68        131.206         131.496
 16          145                   16.31        143.67          143.91
 17          149.66                16.75        148.283         148.493
 18          154.21                15.58        152.788         153.078
 19          163.18                16.51        161.668         161.918
 20          167.58                15.38        166.014         166.254

更快的 Gauss Jacobi 实现不仅明显快于所有其他实现，而且它似乎不像其他方法那样随着数组大小的增加而增加。

在检查正在运行的方法时，似乎快速方法在其第一次迭代时做出了非常好的猜测。

我的猜测是它必须用 np.dot 函数做一些事情，但我不明白为什么这与独立地做每个点积的工作方式不同。

Answer 1

您的第二个实现在每次 k 增量时执行 N 实际 次迭代，因为对 x 的赋值已经覆盖整个向量。它的“优势”因此随着问题的大小而增加。