Участник: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.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-й итерации, что потребует какого-нибудь способа прервать последующие итерации. Для избегания лишних вычисление возможно оптимально оставить эту часть последовательной.