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

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

Основные авторы описания: Анюшева Ирина (614 группа, разделы: 1.1-1.2, 1.7-1.8, 2.7, 3), Исхаков Эльдар (614 группа, разделы 1.3-1.6, 1.9-1.10)

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

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

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

Входные данные алгоритма: система линейных алгебраических уравнений, которую необходимо решить.

Параметры алгоритма: [math]\displaystyle{ m }[/math] - размерность подпространства Крылова, [math]\displaystyle{ \varepsilon }[/math] - параметр останова алгоритма.

Выходные данные: решение системы в виде вектора.

Краткое описание алгоритма:

• Задать параметры алгоритма и начальное приближение;
• Пока не выполнен критерий останова повторить:
  • Построить базис подпространства Крылова;
  • Минимизировать невязку.

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

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

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

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

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

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

В качестве нормы рассматривается евклидова норма.

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

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

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

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

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

Невязку уравнения можно определить как [math]\displaystyle{ 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]\displaystyle{ y - }[/math] вектор размерности [math]\displaystyle{ m }[/math].

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

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

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

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

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

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

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

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

Нижеперечисленные шаги алгоритма и входящие в них операции являются наиболее ресурсоемкими. Они выполняются в приведенном порядке. Каждый из них допускает распараллеливание:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Информационный граф

Метод состоит из двух блоков: ортогонализация и минимизация невязки. При ортогонализации основными операциями являются умножение матрицы на вектор и скалярное умножение векторов, для которых возможен параллелизм, что показано в блоке на темно-синем фоне. В блоке [math]\displaystyle{ min }[/math] решается задача минимизации невязки, которая, строго говоря, является отдельной самостоятельной задачей.

Рис.1. Информационная структура алгоритма GMRES

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

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

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

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

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

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

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

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

• Вычислительная мощность: весьма маленькая и равна [math]\displaystyle{ m^2 }[/math], так как [math]\displaystyle{ m\lt \lt n }[/math].

• Устойчивость: GMRES в этом плане является очень хорошим методом. Теоретически можно получить точное решение не более чем за [math]\displaystyle{ n }[/math] итераций. Однако, очень приемлемые результаты получаются и за [math]\displaystyle{ m\lt \lt n }[/math] итераций. Однако, есть примеры, в которых видно, что первые [math]\displaystyle{ m-1 }[/math] итераций дают невязку, равную некой константе, а на итерации под номером [math]\displaystyle{ m }[/math] получаем точное решение.

Программная реализация алгоритма

Масштабируемость алгоритма и его реализации

Исследование масштабируемости параллельной реализации итерационного метода решения системы линейных алгебраических уравнений GMRES проводилось в тестовой среде суперкомпьютера "Ломоносов". Алгоритм реализован на языке С с использованием библиотеки IntelMPI версии 4.1.0. Сборка производилась при помощи компилятора, разработанным компанией Intel версии 13.1.0. Набор и границы значений изменяемых параметров запуска:

• число процессов [1 : 128] со значениями, равными степеням двойки;
• размер матрицы [128 : 2048] с шагом 128.

Использованные команды для компиляции и запуска:

   mpicc main3.c -o prog

   #!/bin/bash 
   for a in `seq 0 7`; 
   do 
   let b=2**$a 
   sbatch -p test -n $b impi ./prog
   sleep 600s 
   done

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

• минимальная эффективность реализации 0.002%;
• максимальная эффективность реализации 4.92%.

На следующих рисунках приведены графики производительности и эффективности выбранной реализации в зависимости от изменяемых параметров запуска. Размерность подпространства Крылова во всех случаях равна 10.

Рис.2 Изменение производительности выполнения в зависимости от числа процессов и размера входной матрицы

Рис.3 Изменение эффективности выполнения в зависимости от числа процессов и размера входной матрицы

Из рисунков видно, что при увеличении количества процессов, а также при увеличении размерности задачи, производительность параллельной реализации алгоритма уменьшается, что говорит о плохой масштабируемости. Исследуемая реализация алгоритма доступна по ссылке.

Существующие реализации алгоритма

• Реализация алгоритма в библиотеке SciPy [1]
• Реализация алгоритма в библиотеке IML++ [2]
• Библиотека для MATLAB [3]

Литература

1. Y. Saad, M.H. Schultz. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J. Sci. Statist. Comput., 7 (1986), pp. 856–869

2. Е.Е. Тыртышников. Методы численного анализа. – Москва: Академия, 2007

3. М.Ю. Баландин, Э.П. Шурина. Некоторые оценки эффективности параллельных алгоритмов решения СЛАУ на подпространствах Крылова. - Вычислительные технологии, том 3, №1, 1998