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

Основные авторы описания: А.А.Ульянов (разделы 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.9, 2.7), К.Петров (разделы 1.7, 1.8, 1.10, 2.4).

Обобщенный метод минимальных невязок
Последовательный алгоритм
Последовательная сложность	[math] O(nm^2) [/math]
Объём входных данных	[math] (n + 1)^2 [/math]
Объём выходных данных	[math] n [/math]

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

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

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

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

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

• система линейных алгебраических уравнений вида [math] Ax = b [/math], где [math] A [/math] — невырожденная матрица размера [math]n[/math]-на-[math] n [/math];
• [math] x_0 [/math] - начальное приближение.

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

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

Метод GMRES приближает точное решение исходной системы [math] Ax = b [/math] вектором [math] x_m [/math], который минимизирует Евклидову норму невязки [math] r_m = Ax_m - b[/math] на подпространстве Крылова:

[math] K_m = K_m(A,b) = \operatorname{span} \, \{ b, Ab, A^2b, \ldots, A^{m-1}b \} [/math], где [math] \operatorname{span} \, \{ a_1, a_2, \ldots, a_m \} [/math] - линейная оболочка системы векторов [math] \{ a_1, a_2, \ldots, a_m \} [/math].

На каждой итерации метода решение уточняется поправкой, представленной в виде разложения по ортонормированному базису пространства [math] K_m [/math]:

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

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

Для построения ортонормированного базиса [math] v_1, v_2, \ldots, v_m [/math] пространства Крылова [math] K_m [/math] метод использует процесс Арнольди.

Идея процесса Арнольди состоит в том, что для нахождения вектора [math] v_{j+1} [/math] при известных векторах [math] v_1, v_2, \ldots, v_j [/math] необходимо вычислить вектор [math] Av_j [/math] и ортонормировать этот вектор относительно [math] v_1, v_2, \ldots, v_j [/math]:

[math] v_{j+1} = \frac{\hat{v}_{j+1}}{\|\hat{v}_{j+1}\|_2} [/math], где [math] \hat{v}_{j+1} = Av_j - \sum_{i=1}^j h_{ij}v_{i} [/math].

Коэффициенты [math] h_{ij} = (Av_j, v_i) [/math] определяются из условия ортогональности:

[math] A(v_j, v_k) - \sum_{i=1}^j h_{ij}(v_{i}, v_k) = 0 [/math], где [math] \quad k=1,\ldots,j [/math].

Результат процедуры Арнольди совпадает с результатом процесса Грамма-Шмидта, применённого к векторам [math] \{b, Ab, A^2b, \ldots, A^{m-1}b\} [/math].

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

• [math] v_1 [/math], такой что [math] \|v_1\|_2 = 1 [/math];
• матрица [math] A [/math] размером [math] n [/math]-на-[math] n [/math].

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

• [math] v_1, v_2, \ldots, v_m [/math] - ортонормированный базис пространства Крылова [math] K_m [/math].

Алгоритм процесса:

Выполнять для [math] \quad j=1,\ldots,m [/math]:

1. [math] h_{ij} = (Av_j, v_i), \quad i=1,\ldots,j [/math];
2. [math] \hat{v}_{j+1} = Av_j - \sum_{i=1}^j h_{ij}v_{i} [/math];
3. [math] h_{j+1j} = \|\hat{v}_{j+1}\|_2 [/math];
4. [math] v_{j+1} = \frac{\hat{v}_{j+1}}{h_{j+1j}} [/math].

При введении матричных обозначений можно записать:

Поправка решения:

[math] z_m = V_my_m [/math], где [math] y_m \in \mathbb{R}^m [/math] - вектор коэффициентов, [math] V_m = [v_1|v_2|...|v_m] [/math].

Матричное представление процесса Арнольди:

[math] AV_m = V_{m+1}H_{m+1m} = V_mH_m + h_{m+1m}v_{m+1} e_m^T [/math],

где

[math]H_{m+1m} = \begin{bmatrix}h_{11} & h_{12} & h_{13} &... & h_{1m-1} & h_{1m} \\ h_{21} & h_{22} & h_{23} & ... & h_{2m-1} & h_{2m} \\ 0 & h_{32} & h_{33} &...& h_{3m-1} & h_{3m} \\ ... & ... & ... & ... & ... & ... \\ 0 & 0 & 0 & ... & h_{m m-1} & h_{mm} \\0 & 0 & 0 & ... & 0 & h_{m+1 m}\end{bmatrix}[/math]

- верхняя матрица Хессенберга размерности [math] m+1 [/math] на [math] m [/math].

1.2.2 Минимизация функционала невязки

Невязка исходной системы имеет вид:

[math] r_m = b - A(x_0 + z_m) = r_0 - AV_my_m = \beta v_1 - V_{m+1}H_{m+1m}y_m[/math], где [math] \beta = \|r_0\|_2[/math].

C учетом того, что матрица [math] V_{k+1} [/math] составлена из ортонормированных векторов:

[math] r_m = \beta V_{m+1}e_1 - H_{m+1m}y_m[/math], где [math] e_1 = (1,0,0,\ldots,0)^T [/math].

Норма невязки [math] r_m [/math]:

[math] \|r_m\|_2 = \|\beta e_1 - H_{m+1m}y_m\|_2 [/math].

Норма [math] \|r_m\|_2 [/math] достигает минимума при [math] \beta e_1 = H_{m+1m}y_m [/math].

Для нахождения решения такой системы, матрица [math] H_{m+1m} [/math] приводится к верхнему треугольному виду с помощью [math] m [/math] последовательных вращений Гивенса.

1.2.3 Общая схема метода

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

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

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

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

• Вычисление ортонормированного базиса [math] K_m [/math];
• Формирование приближенного решения [math] x_m [/math].

Для вычисления ортонормированного базиса [math] K_m [/math] метод использует процесс Арнольди:

[math] h_{j+1j}v_{j+1} = Av_j - \sum_{i=1}^j h_{ij}v_{j} [/math], где [math] h_{ij} = (Av_j, v_i)[/math], [math] \quad j=1,\ldots,m [/math].

Для нахождения приближённого решения [math] x_m [/math] метод использует формулу:

[math] x_m = x_0 + V_my_m [/math].

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

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

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

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

На практике используется перезапускаемая версия метода - GMRES(m):

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

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

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

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

2 Построение ортонормированного базиса [math] K_m [/math]:

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

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

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

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

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

3 Вычисление [math] x_m [/math]:

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

3.2 Вычислить [math] r_m [/math];

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

4 Перезапуск:

4.1 [math]x_0 = x_m[/math];

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

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

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

Если пренебречь сложностью вычисления [math] y_m [/math], то общую сложность алгоритма GMRES можно разделить на две части:

1. Сложность вычисления ортонормированного базиса пространства [math] K_m [/math]:

a) Для вычисления [math] j [/math]-го вектора базиса [math] K_m, j \lt m [/math] требуется:

[math] NZ + (2j + 1)n [/math] мультипликативных операций, где [math] NZ [/math] - количество ненулевых элементов матрицы [math] A [/math].

b) Вычисление последнего вектора базиса требует:

[math] n(m + 1) [/math] мультипликативных операций.

c) Общая мультипликативная сложность вычисления ортонормированного базиса [math] K_m [/math]:

[math] O(nm^2) [/math].

2. Сложность вычисления приближённого решения [math] x_m [/math]:

a) Вычисление этой формулы требует [math] nm [/math] мультипликативных операций.

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

Рис. 1 input - входные данные алгоритма (квадратная матрица, вектор правой части, начальное приближение, точность решения); [math] A_i [/math] - шаги алгоритма Арнольди, [math] x_i [/math] - нахождение невязки на [math] i [/math]-ой итерации; if - вычисление [math] r_m [/math] и в случае достижения нужной точности - остановка.

Шаги алгоритма Арнольди выполняются строго последовательно, но вычисление [math] i [/math]-ой невязки может происходить параллельно с [math] (i + 1) [/math]-ым шагом алгоритма Арнольди.

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

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

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

При выполнении параллельного алгоритма решения систем линейных алгебраических уравнений необходимо распределить входные данные наиболее рациональным способом. Матрица коэффициентов разбивается на непрерывные горизонтальные полосы. Распределение вектора правой части и начального приближения выполнено по принципу дублирования, т.е. все элементы векторов скопированы на все процессоры. Для рассматриваемой задачи решения СЛАУ такой подход распределения допустим, т.к. векторы содержат столько же элементов, сколько и одна строка или один столбец матрицы коэффициентов.

Подзадачи можно разделить на 2 группы: операции между матрицей и вектором (например, умножение матрицы на вектор) и операции между двумя векторами (например, скалярное умножение). Подзадачи умножения константы на вектор и определение нормы вектора можно относить ко второй группе, так как их трудоемкость одного порядка.

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

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

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

• [math] A [/math] – квадратная матрица размера [math] n [/math]-на-[math] n [/math];
• [math] b [/math] – вектор правой части размера [math] n [/math];
• [math] x_0 [/math] – начальное приближение размера [math] n [/math];
• [math] \epsilon [/math] – погрешность решения.

Для перезапускаемой версии алгоритма дополнительно:

• [math] m [/math] – число итераций.

Объём входных данных: [math] n^2 + 2n + 1 [/math].

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

• [math] x_m [/math] - вектор приближённого решения размера [math] n [/math].

Объём выходных данных: [math] n [/math].

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

Обобщенный метод минимальных невязок предназначен для решения несимметричных систем линейных уравнений. GMRES использует перезазапуски для контроля требований к хранению данных.

Если используется неперезапускаемый метод, GMRES будет сходиться за не более [math] n [/math] шагов. Но это не имеет практического применения, когда [math] n [/math] велико. Кроме того, требования к памяти и вычислениям при отсутствии перезапусков являются очень высокими. Важнейшим элементом для успешного применения [math] GMRES(m)[/math] являются перезапуски; то есть выбор [math] m [/math]. Существуют примеры, для которых метод застаивается и сходимость есть только на [math]m[/math]-м шаге. Для таких систем, любой выбор [math] m [/math] меньше, чем [math] n [/math] - не сходится.

Основным недостатком GMRES является то, что объемы работ и ресурсы памяти, требуемые для итерации, линейно возрастают с увеличением числа итераций. Если не удалось получить чрезвычайно быструю сходимость - затраты будут очень высокими. Основной способ преодолеть это ограничение - это перезапуск итерации. После выбранного числа итераций [math] m [/math] накопленные данные будут очищены и промежуточные результаты используются в качестве входных данных для последующих итераций. Эта процедура повторяется до тех пор, пока сходимость не будет достигнута. Трудность заключается в выборе подходящего значения [math] m [/math]. Если [math] m [/math] слишком мало, [math]GMRES(m)[/math] может быть медленно сходиться, или не сойтись совсем. Значение [math] m [/math], большее, чем это необходимо, включает в себя чрезмерную работу и использует больше памяти. К сожалению, нет определенных правил, регулирующих выбор [math] m [/math]; выбирая , когда делать перезапуск - это вопрос опыта.

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

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

Масштабируемость рассматривалась в работах и . Для исследования масштабируемости была использована параллельная реализация с использованием объекта KSP(решатель СЛАУ) бибилиотеки PETSc(содержащей структуры, типы данных и подпрограммы для параллельной реализации методов на подпространствах Крылова), найденная на GitHub .

Вычисления проводились на вычислительном кластере, состоящем из 16 узлов, со следующей конфигурацией: Altus 1804i Server - 4P Interlagos Node, Quad AMD Opteron 6276, 16C, 2.3 GHz, 128GB, DDR3-1333 ECC, 80GB SSD, MLC, 2.5" HCA, Mellanox ConnectX 2, 1-port QSFP, QDR, memfree, CentOS, Version 5.

На рис. 2 представлен график сильной масштабируемости.

Рис. 2 Зависимость производительности от числа процессоров (запуск производился на 1, 8, 27, 64, 216, 512, 1000 процессорах)

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

Реализации алгоритма GMRES присутствуют во многих программных пакетах:

• LINPACK
• JAMA/C++
• ATLAS
• Intel MLK
• IML++
• LAPACK
• MTL
• SPARSKIT
• WSMP

3 Литература

<references \>