Участник:Lvs/Итерационный метод решения системы линейных алгебраических уравнений GMRES

Эта работа успешно выполнена Преподавателю: в основное пространство, в подстраницу Данное задание было проверено и зачтено. Проверено Evgeny Mortikov и Dan.

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

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

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

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

Дана система линейных алгебраических уравнений [math]Ax = b[/math], где [math]A[/math] - квадратная матрица размера [math]n \times n[/math].

Вычисляется последовательность приближённых решений уравнения [math] x_m [/math], которая сходится к точному решению уравнения [math]x^*[/math].

Подпространство Крылова размерности [math]m[/math], порожденное вектором [math]b[/math] и матрицей [math]A[/math] – линейное пространство [math]K_m(A, v)[/math], являющееся линейной оболочкой векторов [math]b, Ab, A^2b, ..., A^{m-1}b[/math].

Приближённое решение считается по формуле [math] x_m = Q_my_m[/math], где [math]Q_m[/math] - матрица размера [math]n \times m[/math], составленная из векторов ортонормированного базиса пространства [math]K_m(A, v)[/math], [math]y_m \in \mathbb{R}^m[/math].

Итоговым считается решение, минимизирующее норму невязки ([math]||r||_2 = ||b - Ax_m||_2[/math]) с заданной точностью.

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

1. Значение первого вектора выбирается как [math]q_1 = b / ||b||[/math]
2. Оставшиеся векторы [math]q_j, j = 2 .. n[/math] считаются по такому плану:
   1. Вычисляется вектор [math]z[/math] путём умножения матрицы [math]A[/math] на предыдущий вектор [math]q[/math]: [math]z = Aq_{j-1}[/math].
   2. Вычисляются значения [math]h_{ij-1}[/math] для всех [math]i = 1 .. j - 1[/math] как скалярное произведение вектора [math]z[/math] на вектор [math]q_i[/math]: [math]h_{ij-1} = (z, q_i)[/math]. После каждого такого вычисления значение вектора [math]z[/math] понижается на [math]h_{ij-1} q_i[/math]: [math]z = z - h_{ij-1} q_i[/math]
   3. Затем вычисляется значение [math]h_{jj-1}[/math] как норма вектора [math]z[/math]: [math]h_{jj-1} = ||z||[/math]. Алгоритм останавливается, если эта норма равна 0
   4. Если полученное значение нормы не равно 0, то [math]q_j = z /||z||[/math] (что эквивалентно [math]q_j = z /h_{jj-1}[/math]

Сложность алгоритма [math]O(n^2)[/math].

Получаемые в процессе работы алгоритма векторы [math]q_1, q_2, ..., q_m[/math] образуют искомый ортонормированный базис Крыловского пространства [math]K_m(A, b)[/math].

Значения [math]h_{i,j}[/math] образуют верхнюю матрицу Хессенберга (обозначим как [math]H_m[/math], размерность [math]m+1 \times m[/math]).

[math]H_m = \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]r_m = H_my_m - ||b - Ax_0||e_1[/math], где [math]x_0[/math] - вектор начального приближения. Минимизация второй нормы вектора невязки представляет из себя решение линейной задачи минимальных квадратов размера [math]n[/math].

Приближённое решение вычисляется по формуле [math] x_m = Q_my_m[/math] после нахождения соответствующего значения минимизирующего вторую невязку вектора. Сложность вычислений [math]O(nm)[/math]

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

Самая трудоёмкая и самая часто вызываемая операция в алгоритме - умножение матрицы на вектор. В первую очередь это связано с вычислением матрицы Хессенберга и ортонормированного базиса подпространства Крылова. Число умножений матрицы на вектор: [math]n-1[/math] в алгоритме Арнольди, K в выбранном алгоритме решения задачи минимальных квадратов и 1 в подсчёте результата. Также в алгоритм входят [math]n*(n-1) / 2[/math] операций скалярного произведения.

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

1. Алгоритм Арнольди
2. Решение линейной задачи минимальных квадратов
3. Вычисление решения

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

    Вычислить ортонормированный базис размера m
    Вычислить минимальное приближение (минимальную норму вектора невязки)
    if (приближение < заданная точность)
        Вычислить результат
        return;
    else
        Повторить для m = m+1

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

Сложность вычислений собственных векторов на каждом шаге - [math]O(n^2)[/math].

Сложность решения задачи минимальных квадратов [math]O(m^2)[/math], но при последовательном выполнении алгоритма можно понизить сложность до [math]O(m)[/math] путём использования подсчитанных на прошлых итерациях значений.

Сложность вычисления решения [math]O(nm)[/math], поэтому его целесообразно вычислять только в конце.

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

Основной ресурс параллелизма - умножение матрицы на вектор, поскольку является самым ресурсоёмким, но при этом умножать отдельные строки матрицы можно независимо друг от друга. Для этой операции информационный граф представлен на рисунке.

Рис.1. Информационный граф умножения матрицы на вектор

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

1. Параллельное умножение матрицы на вектор выполняется в [math] n [/math] ярусов умножений и сложений.
2. Операция вычисления скалярного произведения
3. Шаги алгоритма Арнольди выполняются строго последовательно, однако последующее нахождение невязки можно выполнять параллельно с запуском следующего шага алгоритма.
4. Сложность решения задачи минимальных квадратов зависит от используемого алгоритма. В отличие от последовательного метода данные, полученные на прошлых шагах здесь использовать нельзя.
5. Вычисление результата выполняется один раз, его можно выполнить параллельно в [math] n [/math] ярусов умножений и сложений.

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

Входные данные: Матрица [math]A[/math], вектор [math]b[/math], точность [math]t[/math]. Объём: [math]n^2 + n + 1[/math]

Выходные данные: Вектор [math]x_m[/math]. Объём: [math]n[/math]

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

Теоретически на каждой итерации алгоритма сложность вычислений составляет [math]n + X[/math] операций умножения, где [math]X[/math] - сложность алгоритма для решения задачи минимальных квадратов. Таким образом в худшем случае сложность будет примерно равна [math]n^2 + X*n + n[/math] операций умножения.

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

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

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

Для исследования масштабируемости алгоритма была использована параллельная реализация, найденная на просторах GitHub. Эта реализация использует многопроцессорное выполнение с общей памятью и вычисляет параллельно скалярное произведение и умножение матрицы на вектор. Программа была запущена на суперкомпьютере BlueGene на разреженной матрице размера 90448х90448 с 1921955 ненулевыми элементами, расположенными около главной диагонали. Параметрами запуска являются 100 итераций алгоритма и 10 перезапусков алгоритм (повторный запуск алгоритма с использованием решения, полученного на предыдущем запуске, в качестве начального приближения). График зависимости времени выполнения от числа процессоров представлен на рисунке 2.

Рис 2. График зависимости времени выполнения от числа процессоров

График эффективности для полученного времени выполнения указан на рисунке 3.

Рис 3. График зависимости эффективности реализации от числа процессоров

Параметры компиляции: mpixlc_r -qsmp=omp -O3

Параметры запуска: gmres mat.mtx void 100 0 10

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

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

Реализации алгоритма GMRES присутствуют во многих программных продуктах как части стандартных библиотек, например Matlab или DEFORM3D. Также существую отдельные библиотеки, например MGMRES для языка С++.

3 Литература

[1] Тыртышников Е. Е. Методы численного анализа — М., Академия, 2007.

[2] М. Дана, А. Г. Зыков, Х. Д. Икрамов Метод минимальных невязок для специального класса линейных систем с нормальными матрицами коэффициентов - Ж. вычисл. матем. и матем. физ., 2005, том 45, номер 11