Участник: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.3 Вычислительное ядро алгоритма

Алгоритм Арнольди:

q[1] = b/||b||
for j = 1, ..., n-1 do
    z = Aq[j]; //очередной вектор
    for i = 1,..., j do
        h[i][j] = (z, q[i]);
        z = z − h[i][j]q[i];
    end i
    h[j+1][j] = ||z||;
    if (h[j+1][j] == 0) break;
    q[j+1] = z/h[j+1][j];
end j

Сложность алгоритма [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_k = \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.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.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.7 Существующие реализации алгоритма

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