Участник:Анюшева Ирина/Итерационный метод решения системы линейных алгебраических уравнений GMRES (обобщенный метод минимальных невязок)

GMRES (обобщенный метод минимальных невязок)
Последовательный алгоритм
Последовательная сложность	[math]O(m^2n^2)[/math]
Объём входных данных	1
Объём выходных данных	1
Параллельный алгоритм
Высота ярусно-параллельной формы	1
Ширина ярусно-параллельной формы	1

Основные авторы описания: Анюшева Ирина (614 группа), Исхаков Эльдар (614 группа)

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

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

Обобщенный метод минимальных невязок (GMRES) для решения системы линейных алгебраических уравнений аппроксимирует решение с помощью вектора в подпространстве Крылова с минимальным остатком. Метод GMRES разработан Ю.Саадом и Мартин Х. Шульцом в 1986 году. Метод по праву считается одним из самых эффективных численных методов решения несимметричных систем. Он гарантирует неувеличение нормы невязки в ходе итерационного процесса, который основан на построении базиса в соответствующей системе подпространстве Крылова [math]K_j(A, r_0)[/math]. Затем решение уточняется некоторой добавкой, представленной в виде разложения по этому базису.

Метод обобщенных минимальных невязок популярен из-за наличия ряда преимуществ: он ошибкоустойчив, допускает эффективное распараллеливание, не требует нахождения параметра релаксации, обладает суперлинейной скоростью сходимости. Однако в чистом виде для больших систем метод не используется из-за чрезвычайно больших затрат памяти. Зато широкое распространение получила перезапускаемая версия метода, подоразумевающая перезапуск метода каждые m итераций. Уменьшение размерности пространства может привести к стагнации метода – процессу, при котором с каждой следующей итерацией норма невязки уменьшается незначительно, либо не уменьшается вовсе.

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

Исходные данные: система линейных алгебраических уравнений [math]Ax=b[/math], где:

[math]A=(a_{ij})[/math] [math]-[/math] вещественная матрица размера [math]n \times n[/math];

[math]b[/math] и [math]x[/math] вектора из [math]n[/math] элементов;

[math]x^*[/math] [math]-[/math] точное решение системы.

Вычисляемые данные: [math]x^{(s)} - [/math]приближенное решение системы.

1.2.1 Ортогонализация Арнольди

Для построения базиса [math]V[/math] в подпространстве Крылова в методе GMRES применяется ортогонализация Арнольди.

Формулы метода:

[math] \begin{align} (h_{i+1,i}v_{i+1}) & = Av_i - \sum_{j = 1}^{i} h_{j,i}v_j, \\ v_1 &= \frac{r_0}{\parallel r_0 \parallel}, \end{align} [/math]

В матрично-векторных обозначениях соотношение может быть записано как

[math] \begin{align} AV_i &= V_i H_i + h_{i+1,i}v_{i+1} e_i^T &=V_{i+1}\overline{H}_i, \end{align} [/math]

[math] \begin{align} V_i &= [v_1|v_2|...|v_i], \overline{H}_i - \end{align} [/math] верхняя матрица Хессенберга размерности [math](i+1)\times i[/math]

1.2.2 Минимизация невязки

Невязку уравнения можно определить как [math] J(y) = \parallel b- A x\parallel = \parallel b- (Ax_0+V_m y)\parallel =\parallel \beta e_1 - \overline{H}_m y \parallel[/math], где [math]y -[/math] вектор размерности [math]m[/math].

Вектор [math]x_m[/math] может быть получен как [math]x_m = x_0 + V_m y_m[/math], вектор [math]y_m[/math] может найден как решение линейной задачи наименьших квадратов размера [math](m+1)\times m[/math], где [math]m\lt \lt n[/math].

Для решения задачи минимизации приведем матрицу [math]\overline{H}_m[/math] к верхнему треугольному виду с помощью вращений Гивенса. На шаге [math]i[/math] матрица вращения [math]\Omega_i[/math] имеет размеры [math](m+1)\times (m+1)[/math] и имеет следующий вид:

[math] \Omega_i= \begin{pmatrix} 1 & &&&&&& \\ &\ddots &&&&&&\\ &&1&&&&& \\ &&&c_i&s_i &&& \\ &&&-s_i& c_i &&&\\ &&&&&1&&& \\ &&&&&&\ddots&\\ &&&&&&&1 \\ \end{pmatrix}[/math]

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

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

• умножение исходной матрицы на вектор при вычислении невязок;

• вычисление скалярных произведений векторов при построении ортонормированного базиса подпространства Крылова.

1.4 Макроструктура алгоритма

Наиболее ресурсоемкие операции при выполнении алгоритма:

• Вычисление нормы невязки для задачи "начальных" условий;

• Вычисление элементов Хессенберговой матрицы - матрично-векторные умножения и скалярные произведения векторов;

• Вычисление произведений разреженных матриц при минимизации невязки.

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

1.5.1 Шаги алгоритма

1. Выбрать произвольное [math]x_0[/math]; вычислить [math]r_0=b-Ax_0[/math] и [math]v_1=\frac{r_0}{\parallel r_0 \parallel}[/math]

2. Задать матрицу [math]\overline{H}_m[/math] размерности [math](m+1)\times m[/math] и всем ее элементам присвоить нулевые значения

3. Для [math]j=1,\dots,m[/math] выполнять:

• для [math]i=1,\dots,j[/math] выполнять: [math]h_{i,j}=(Av_j,v_i)[/math]
• [math]\omega_j=Av_j-\sum_{i=1}^j h_{i,j}v_i[/math]
• [math]h_{j+1,j}={\parallel \omega_j \parallel}[/math]
• [math]v_{j+1}=\frac{\omega_j}{h_{j+1,j}}[/math]

4. Вычислить [math]y_m=\arg\min_{y} \parallel \beta e_1 - \overline{H}_my \parallel,\ y \in \mathbb{R}^m[/math] и [math]x_m=x_0+V_my_m[/math]

5. Вычислить [math]r_m=b-Ax_m[/math]. Если достигнутая точность удовлетворительна, остановиться.

Иначе принять [math]x_0=x_m[/math] и [math]v_1=\frac{r_m}{\parallel r_m \parallel}[/math] и вернуться к пункту 2.

1.5.2 Решение задачи минимизации невязки

Для нахождения [math]y_m=\arg\min_{y} \parallel \beta e_1 - \overline{H}_my \parallel[/math] необходимо выполнить следующие шаги:

1. Найти [math]m[/math] матриц вращений Гивенса [math]\Omega_i(c_i,s_i),\ i=1,\dots,m[/math] размера [math](m+1)\times(m+1)[/math], где [math]s_i=\frac{h_{i+1,i}}{\sqrt{h_{i,i}^2+h_{i+1,i}^2}}[/math], [math]c_i=\frac{h_{i,i}}{\sqrt{h_{i,i}^2+h_{i+1,i}^2}}[/math]

2. Вычислить [math]\overline{H}_m^{(m)}=\Omega_1\dots\Omega_m\overline{H}_m[/math] и [math]\overline{g}_m=\Omega_1\dots\Omega_m\beta e_1[/math], удалить последнюю строку и последний элемент соответственно и решить СЛАУ с полученной верхнетреугольной матрицей и правой частью. Полученное решение [math]y_m[/math] минимизирует исходный функционал.

1.6 Последовательная сложность алгоритма

Для выполнения одной итерации алгоритма необходимо произвести:

• [math]3m+1[/math] делений;

• [math]2m+1[/math] вычислений квадратного корня;

• [math]n^2+\frac{3m+5}{2}n[/math] умножений;

• [math]O(\frac{m^2n^2}{2})[/math] сложений (вычитаний).

1.8 Ресурс параллелизма алгоритма

Для нахождения матрично-векторных произведений при вычислении невязок матрица СЛАУ должна быть подвергнута декомпозиции на строчные блоки, размеры которых определяются требованиями равномерности загрузки. Соответствующим образом подвергается декомпозиции и вектор правой части. Для наиболее эффективного вычисления скалярных произведений и линейных комбинаций векторов базиса каждый из векторов [math]v_j[/math] должен быть разбит на блоки равного размера по числу процессов; каждый из процессов при этом хранит соответствующие блоки всех векторов. Обработка вектора решения дублируется на всех процессах.

1.9 Входные и выходные данные алгоритма

Входные данные:

• матрица [math]A[/math] размера [math]n\times n[/math];
• вектор-столбец [math]b[/math] размера [math]n[/math];
• вектор-столбец [math]x_0[/math] размера [math]n[/math] - начальное приближение;
• натуральное число [math]m\lt \lt n[/math], определяющее размерность подпространства Крылова.

Выходные данные:

• вектор-столбец [math]x[/math] размера [math]n[/math] - решение системы.