Что случилось с моей реализацией GMRES?

Question

Я пытаюсь реализовать GMRes в Jupyter Notebook, которая (в случае, если вы не знаете):

Это мой код:

import numpy as np 

def GMRes(A, b, x0, e, nmax_iter, restart=None): 
    r = b - np.asarray(np.dot(A, x0)).reshape(-1) 

    x = [] 
    q = [0] * (nmax_iter) 

    x.append(r) 

    q[0] = r/np.linalg.norm(r) 

    h = np.zeros((nmax_iter + 1, nmax_iter)) 

    for k in range(nmax_iter): 
     y = np.asarray(np.dot(A, q[k])).reshape(-1) 

     for j in range(k): 
      h[j, k] = np.dot(q[j], y) 
      y = y - h[j, k] * q[j] 
     h[k + 1, k] = np.linalg.norm(y) 
     if (h[k + 1, k] != 0 and k != nmax_iter - 1): 
      q[k + 1] = y/h[k + 1, k] 

     b = np.zeros(nmax_iter + 1) 
     b[0] = np.linalg.norm(r) 

     result = np.linalg.lstsq(h, b)[0] 

     x.append(np.dot(np.asarray(q).transpose(), result) + x0) 

    return x 

По мне это должна быть правильной, но когда я исполню:

A = np.matrix('1 1; 3 -4') 
b = np.array([3, 2]) 
x0 = np.array([1, 2]) 

e = 0 
nmax_iter = 5 

x = GMRes(A, b, x0, e, nmax_iter) 

print(x) 

Примечание: На данный момент e ничего не делает.

я получаю это:

[array([0, 7]), array([ 1., 2.]), array([ 1.35945946, 0.56216216]), array([ 1.73194463, 0.80759216]), array([ 2.01712479, 0.96133459]), array([ 2.01621042, 0.95180204])] 

x[k] должны приближаться к (32/7, -11/7), так как это результат, но вместо этого он приближается к (2, 1), что я делаю не так?

Answer 1

Я думаю, что алгоритм дает правильный результат.

Вы пытаетесь решить Ах = Ь, где:

Если попытаться найти решение вручную, операция матрицы вы пытаетесь решить эквивалентна системе, которая может быть решена используя замену.

Если вы пытаетесь решить, вы увидите, что решение:

что то же самое решение, ваш алгоритм дает.

Вы можете дважды проверить это с помощью GMRES реализации уже в SciPy:

import scipy.sparse.linalg as spla 
import numpy as np 

A = np.matrix('1 1; 3 -4') 
b = np.array([3, 2]) 
x0 = np.array([1, 2]) 
spla.gmres(A,b,x0) 

Какие выходы

array([ 2., 1.]) 

Answer 2

Обратите внимание, что этот алгоритм сходящуюся к правильному результату, но делает это слишком медленно , Максимальное число итераций GMRES для сходимости никогда не должно превышать размерность матрицы A. Если размерность матрицы A равна n, то вектор (n + 1) 'th Arnoldi должен быть равен нулю, например. мы должны иметь возможность полностью охватить пространство Крылова с векторами n Арнольди. Я бы просто применить следующий патч, и все должно работать, как они должны: вектор последовательность

- for k in range(nmax_iter): 
+ for k in range(min(nmax_iter, A.shape[0])): 
     y = np.asarray(np.dot(A, q[k])).reshape(-1) 

-  for j in range(k): 
+  for j in range(k + 1): 

Раствор должен затем быть просто:

[array([ 1.  , 0.35294118]), array([ 2., 1.])] 

например мы сходились в двух итерациях, которые мы ожидаем, так как n = 2.