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

1 ЧАСТЬ. Свойства и структура алгоритмов

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

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

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

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

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

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

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

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

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

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

На каждом шаге алгоритма получаются векторы $q_1, q_2, ..., q_m$, являющиеся ортонормированным базисом Крыловского пространства $K_m(A, b)$.

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

$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}$.

Вектор невязки задан уравнением $r_m = H_my_m - ||b - Ax_0||e_1$, где $x_0$ - вектор начального приближения. Минимизация второй нормы вектора невязки представляет из себя решение линейной задачи минимальных квадратов размера $n$.

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

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

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

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

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

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

Сложность вычислений собственных векторов на каждом шаге - $O(n^2)$, но может уменьшиться до $O(n)$ при более разреженной матрице.

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

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

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

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

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

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

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

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

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

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