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

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

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

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

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

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

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

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

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

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

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

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

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]\displaystyle{ O(n^2) }[/math].

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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