为什么GaussJacobi方法的特定numpy实现会显著减少迭代次数?

2024-09-28 17:02:03 发布

您现在位置:Python中文网/ 问答频道 /正文

在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方法。这两种方法的实现方式与我最初的慢高斯-雅可比方法类似

我对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函数有关,但我不明白为什么这与独立处理每个点积不同


Tags: 方法defnp次数dotdiagsolutionshape