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

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

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

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

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 .. m[/math] считаются по такому плану:
   1. Вычисляется вектор [math]z[/math] путём умножения матрицы [math]A[/math] на предыдущий вектор [math]q[/math]: [math]z = Aq_{j-1}[/math].
   2. Вычисляются значения [math]h_{ij}[/math] для всех [math]i = 1 .. j - 1[/math] как скалярное произведение вектора [math]z[/math] на вектор [math]q_i[/math]: [math]h_{ij} = (z, q_i)[/math]. После каждого такого вычисления значение вектора [math]z[/math] понижается на [math]h_{ij} q_i[/math]: [math]z = z - h_{ij} q_i[/math]
   3. Затем вычисляется значение [math]h_{jj-1}[/math] как норма вектора [math]z[/math]: [math]h_{jj-1} = ||z||[/math]. Алгоритм останавливается, если эта норма равна -
   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 Вычислительное ядро алгоритма

Самая трудоёмкая и самая часто вызываемая операция в алгоритме - умножение матрицы на вектор. В первую очередь это связано с вычислением матрицы Хессенберга.

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

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

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

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

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

Сложность вычислений собственных векторов на каждом шаге - [math]O(n^2)[/math], но может уменьшиться до [math]O(n)[/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. Вычисление результата выполняется один раз, его можно выполнить параллельно в [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. Параметрами запуска являются 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