При реализации алгоритма Гаусса Якоби в python я обнаружил, что две разные реализации требуют существенно различного числа итераций, чтобы сходиться.
Первая реализация - это то, что я изначально придумал
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
Вторая реализация основана на этой статье.
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 итераций.
Что делает вторую реализацию намного быстрее, чем первая?
РЕДАКТИРОВАТЬ:
Я реализовал два других метода, метод Гаусса-Зейделя и метод 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, но я не понять, почему это будет работать по-другому, чем делать каждый точечный продукт независимо.