Участник:Alexander34396/Обобщенный метод минимальных невязок

Свойства и структура алгоритма

Общее описание алгоритма

Обобщённый метод минимальных невязок (англ. Generalized minimal residual method, GMRES) - итерационный метод численного решения системы линейных алгебраических уравнений с невырожденной матрицей. Метод основан на минимизации квадратичного функционала невязки на подпространствах Крылова. Разработан Юсефом Саадом и Мартином Шульцем в 1986 году как обобщение метода MINRES на случай систем с несимметричными матрицами.

Математическое описание алгоритма

Исходные данные:

• система линейных алгебраических уравнений вида [math]\displaystyle{ Ax = b }[/math], где [math]\displaystyle{ A }[/math] — невырожденная матрица размера [math]\displaystyle{ n }[/math]-на-[math]\displaystyle{ n }[/math];
• подпространство Крылова размерности [math]\displaystyle{ m, m \lt = n }[/math], порождённое вектором [math]\displaystyle{ b }[/math] и матрицей [math]\displaystyle{ A }[/math]:

[math]\displaystyle{ K_m = K_m(A,b) = \operatorname{span} \, \{ b, Ab, A^2b, \ldots, A^{m-1}b \}. \, }[/math]

Вычисляемые данные:

• [math]\displaystyle{ x_m \in K_m }[/math] - приближённое решение исходной системы.

Метод GMRES приближает точное решение исходной системы [math]\displaystyle{ Ax = b }[/math] вектором [math]\displaystyle{ x_m \in K_m }[/math], минимизирующим Евклидову норму невязки [math]\displaystyle{ r_m = Ax_m-b }[/math].

Для решения исходной системы GMRES, используя [math]\displaystyle{ l_2 }[/math]-ортонормальный базис пространства [math]\displaystyle{ K_m }[/math], выполняет поиск приближённого решения [math]\displaystyle{ x_m }[/math] в виде:

[math]\displaystyle{ x_m = x_0 + z_m }[/math],

где [math]\displaystyle{ x_0 }[/math] - некоторое начальное приближение, [math]\displaystyle{ z_m \in K_m }[/math] - поправка решения.

Для построения ортонормального базиса [math]\displaystyle{ K_m }[/math] метод использует ортогонализацию Арнольди. При введении для базиса [math]\displaystyle{ K_m }[/math] матричного обозначения [math]\displaystyle{ V_m }[/math] можно записать:

[math]\displaystyle{ z_m = V_my_m }[/math],

где [math]\displaystyle{ y_m \in \mathbb{R}^m }[/math] - вектор коэффициентов.

В общем виде алгоритм метода GMRES может быть записан следующим образом:

1. найти ортонормальный базис [math]\displaystyle{ V_m }[/math] подпространства [math]\displaystyle{ K_m }[/math] с помощью ортогонализации Арнольди;
2. найти [math]\displaystyle{ y_m }[/math], минимизирующий [math]\displaystyle{ \|r_m\|_2 }[/math];
3. посчитать [math]\displaystyle{ x_m = x_0 + V_my_m }[/math];
4. посчитать [math]\displaystyle{ r_m }[/math] и остановиться,если требуемая точность была достигнута, иначе повторить.

Вычислительное ядро алгоритма

Вычислительное ядро последовательной версии метода GMRES состоит из двух частей:

• Вычисление ортонормального базиса [math]\displaystyle{ K_m }[/math] с помощью ортогонализации Арнольди:

на одной итерации алгоритма [math]\displaystyle{ m(m + 1) n + mNZ }[/math] мультипликативных операций, где NZ - количество ненулевых элементов матрицы [math]\displaystyle{ A }[/math];

• Формирование приближенного решения [math]\displaystyle{ x_0 + V_mY_m }[/math]:

[math]\displaystyle{ nm }[/math] мультипликативных операций.

Схема реализации последовательного алгоритма

Для решения исходной системы методом GMRES можно воспользоваться следующим алгоритмом:

1 Подготовка перед итерационным процессом:

1.1 Выбрать начальное приближение [math]\displaystyle{ x_0 }[/math];

1.2 Посчитать невязку [math]\displaystyle{ r_0 = b - Ax_0 }[/math];

1.3 Вычислить [math]\displaystyle{ v_1 = \frac{r_0}{\|r_0\|_2} }[/math].

2 Построение ортонормального базиса [math]\displaystyle{ K_m }[/math]:

Для всех [math]\displaystyle{ j }[/math] от 1 до m по нарастанию выполнять:

2.1 [math]\displaystyle{ h_{ij} := (Av_j, v_i), \quad i=1,\ldots,j }[/math];

2.2 [math]\displaystyle{ \hat{v}_{j+1} := Av_j - \sum_{i=1}^j h_{ij}v_{i} }[/math];

2.3 [math]\displaystyle{ h_{j+1j} = \|\hat{v}_{j+1}\|_2 }[/math];

2.4 [math]\displaystyle{ v_{j+1} = \frac{\hat{v}_{j+1}}{h_{j+1j}} }[/math].

3 Вычисление приближённого решения [math]\displaystyle{ x_m }[/math]:

3.1 [math]\displaystyle{ x_m = x_0 + V_my_m }[/math], где [math]\displaystyle{ y_m }[/math] минимизирует [math]\displaystyle{ \|r_0 - AV_my_m\|_2 }[/math];

3.2 Вычислить [math]\displaystyle{ r_m }[/math];

3.3 Если требуемая точность достигнута, остановиться.

4 Рестарт:

4.1 [math]\displaystyle{ x_0 = x_m }[/math];

4.2 [math]\displaystyle{ v_1 = \frac{r_m}{\|r_m\|_2} }[/math];

4.3 Перейти к шагу 2.