Участник:Макарова Алёна/Итерационный метод решения системы линейных алгебраических уравнений GMRES (обобщенный метод минимальных невязок)

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

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

Метод минимальных невязок представляет собой итерационный метод для численного решения несимметричной системы линейных уравнений. Метод приближает решение вектором в подпространстве Крылова с минимальной невязкой. Чтобы найти этот вектор используется итерация Арнольди.

Впервые метод минимальных невязок был описан в книге Марчука и Кузнецова "Итерационные методы и квадратичные функционалы" и представлял собой геометрическую реализацию метода минимальных невязок.

Алгебраическая реализация метода минимальных невязок была предложена Саадом и Шульцем в 1986 году . С их подачи за методом закрепилась аббревиатура GMRES (обобщенный метод минимальных невязок).

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

Для решения системы [math]Ax = b[/math] с невырожденной матрицей [math]A[/math], выбирается начальное приближение [math]x_0[/math], затем решается редуцированная система [math]Au = r_0[/math], где [math]r_0 = b - Ax_0, x = x_0 + u[/math]. Подпространства Крылова строятся для редуцированной системы (в предположении, что [math]r_0 \neq 0[/math]).

Пусть [math]x_i = x_0 + y[/math], где [math] y \in \mathcal{K}_i [/math]. Невязка имеет вид [math]r_i = r_0 - Ay[/math], а её длина минимальна, в силу теоремы Пифагора, в том и только том случае, когда

[math] r_i \perp A\mathcal{K}_i [/math].

Таким образом, для реализации [math]i[/math]-го шага требуется опустить перпендикуляр из вектора [math]r_0[/math] на подпространство [math]A\mathcal{K}_i[/math]. Проще всего это сделать, если в данном подпространстве уже найден ортогональный базис.

1.2.1 Геометрическая реализация

Геометрическая реализация метода минимальных невязок заключается в построении последовательности векторов [math]q_1, q_2,...[/math] таким образом, что [math]q_1, q_2,.., q_i[/math] дают базис в подпространстве Крылова [math]\mathcal{K}_i[/math] и при этом вектора [math]p_1 = Aq_1, ..., p_i = Aq_i[/math] образуют ортогональный базис в подпространстве [math]A\mathcal{K}_i[/math]. Вектор [math]q_{i+1} \notin \mathcal{K}_i[/math] должен обладать следующими свойствами:

[math] q_{i+1} \in \mathcal{K}_{i+1}, p_{i+1} = Aq_{i+1} \perp A\mathcal{K}_i [/math].

Его можно получить с помощью процесса ортогонализации, применённого к вектору [math]p = Aq[/math], где [math]q = Aq_i[/math]. В геометрической реализации необходимо хранить две системы векторов: [math]q_1, ..., q_i[/math] и [math]p_1, ..., p_i[/math].

1.2.2 Алгебраическая реализация

Алгебраическая реализация метода минимальных невязок использует лишь одну систему векторов, образующих ортогональные базисы в подпространствах [math]\mathcal{K}_i[/math].

[math]q_1 = r_0/||r_0||_2[/math]. Чтобы получить ортонормированный базис [math]q_1, ..., q_{i+1}[/math] в [math]\mathcal{K}_{i+1} = \mathcal{K}_{i+1}(r_0, A)[/math], проводим ортогонализацию вектора [math]Aq_i[/math] к векторам [math]q_1, ..., q_i[/math]. Если [math]Q_i \equiv [q_1, ..., q_i] \in \mathbb{C}^{n*i}[/math], то получим

[math] AQ_i = Q_{i+1}\hat{H_i}, \hat{H_i}= \begin{bmatrix} H_i \\ 0 ... 0 & h_{i+1 i} \end{bmatrix}, [/math]

где [math]H_i[/math] - верхняя хессенбергова матрица порядка [math]i[/math].

Далее рассмотрим [math]QR[/math]-разложение прямоугольной матрицы [math]\hat{H_i}[/math]:

[math] \hat{H_i} = U_iR_i, R_i \in \mathbb{C}^{i*i} [/math].

Тогда минимум невязки [math]||r_0 - AQ_iy||_2[/math] по всем [math]y[/math] будет достигаться в том случае, если [math]y = y_i[/math] удовлетворяет уравнению

[math] R_iy_i = z_i \equiv ||r_0||_2U^T_ie_1, e_1 = \begin{bmatrix} 1 0 ... 0 \end{bmatrix}^T, [/math].

Матрица [math]R_i[/math] невырожденная, следовательно,

[math] x_i = x_0 + Q_iy_i = x_0 + Q_iR_i^{-1}z_i [/math].

1.2.3 Сходимость метода

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

Но об этом не верно говорить в общем. Согласно теореме (Greenbaum, Pták and Strakoš), для каждой монотонной последовательности [math] a_1,..., a_{m-1}, a_m = 0 [/math], найдётся матрица [math]A[/math] такая что [math]||r_n|| = a_n[/math] для всех [math]n[/math], где [math]r_n[/math] невязка, определенная ранее.

На практике для матрицы порядка [math]n[/math] число шагов может оказаться равным [math]n[/math], а норма невязки может оставаться равной норме начальной невязки на всех шагах, кроме последнего.

Пример

[math] \begin{bmatrix} 0 1 & ... \\ 0 0 1 & ... \\...\\... & 01 \\1 & ... \end{bmatrix} \begin{bmatrix} 1 \\ 0 \\...\\ 0 \\ 0 \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\...\\ 0 \\ 1 \end{bmatrix} [/math]

Если [math] x_0 = 0 [/math], то [math] r_0 = b[/math]. Подпространства Крылова получаются такие:

[math] \mathcal{K}_i = span \{ e_n, e_{n-1}, ..., e_{n-i+1}\}, i \le i \le n [/math]

Точное решение системы имеет вид [math]x = e_1[/math], а на итерациях получаются следующие векторы:

[math] x_0 = x_1 = ... = x_{n-1} = 0, x_n = e_1 [/math].

Для получения оценок требуются дополнительные предположения относительно матрицы коэффициентов.

1.2.3.1 Сходимость для положительно определенной матрицы

Если матрица [math]A[/math] положительно определена, то

[math] \|r_n\| \leq \left( 1-\frac{\lambda_{\mathrm{min}}^2(1/2(A^T + A))}{ \lambda_{\mathrm{max}}(A^T A)} \right)^{n/2} \|r_0\|, [/math],

где [math]\lambda_{\mathrm{min}}(M)[/math] и [math]\lambda_{\mathrm{max}}(M)[/math] обозначают наибольшее и наименьшее собственное значение матрицы [math]M[/math], соответственно.

1.2.3.2 Сходимость для симметричной матрицы

Если матрица [math]A[/math] является симметрично и положительно определена, тогда

[math] \|r_n\| \leq \left( \frac{\kappa_2(A)^2-1}{\kappa_2(A)^2} \right)^{n/2} \|r_0\| [/math],

где [math]\kappa_2(A)[/math] число обусловленности матрицы [math]A[/math] в Евклидовой норме.

Если не требовать положительной определённости

[math] \|r_n\| \le \inf_{p \in P_n} \|p(A)\| \le \kappa_2(V) \inf_{p \in P_n} \max_{\lambda \in \sigma(A)} |p(\lambda)| \|r_0\|, \, [/math]

где [math]P_n[/math] обозначает множество полиномов степени не выше [math]n[/math] и [math]p(0) = 1[/math], матрица [math]V[/math] появляется из спектрального разложения матрицы [math]A[/math] и [math]\sigma(A)[/math] - спектр матрицы [math]A[/math]. Грубо говоря, это означает, что быстрая сходимость имеет место, когда собственные значения матрицы [math]A[/math] группируются от происхождения и [math]A[/math] не сильно отклоняется от нормальной матрицы (are clustered away from the origin and A is not too far from).

Все эти неравенства связаны с невязкой вместо фактической ошибки, т.е. расстоянием между текущей итерацией и точным решением.

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

Самой трудоёмкой процедурой в методе является [math]QR[/math]-разложение [math]\hat{H_i}[/math]. В общем случае требуется [math]\mathcal{O}(i^3)[/math] арифметических операций, однако то, что [math]\hat{H_i}[/math] является хессенберговой матрицей, позволяет сократить требуемый объём до [math]\mathcal{O}(i^2)[/math] действий.

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

Алгоритм начинается с выбора начального приближения [math]x_0[/math]. Далее для [math]i = 1,2 ..[/math]

1. Вычисляем матрицу [math]\hat{H_i}[/math].

2. Находим её [math]QR[/math]-разложение: [math]\hat{H_i} = U_iR_i[/math].

3. Решаем систему [math]R_iy_i = U^T_ie_1[/math].

4. Находим приближение [math]x_i = x_0 + Q_iy_i[/math].

Алгоритм завершается, когда норма вектора невязки становится достаточно малой. При большом числе итераций вычислительные затраты и требуемая для хранения векторов [math]p_i[/math] память могу оказаться чрезмерно большими. В таких случаях применяют алгоритм с перезапусками (restart): задают максимальную допустимую размерность [math]K[/math] пространств Крылова [math]\mathcal{K}_i[/math], и если точность не достигнута, берут полученное приближение в качестве нового начального приближения [math]x^0_{new} = x_K[/math].

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

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

Будем предполагать, что матрица решаемой СЛАУ является разреженной. Пусть несимметричная матрица [math]A[/math] имеет симметричный портрет и число ненулевых элементов в нижнем(верхнем) треугольнике равно [math]k[/math].

При [math] m \lt \lt n [/math] можно говорить, что один GMRES-цикл имеет сложность [math]\mathcal{O}(nm^2)[/math].

Требования по памяти с учетом [math]1 \lt m \lt \lt n[/math] и [math]k \lt \lt n[/math] составляют [math]\mathcal{O}(mn)[/math]..

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

В вычислительной схеме присутствуют два интенсивных в вычислительном плане этапа:

1. Умножение матрицы СЛАУ на вектор (для нахождения невязок).

Для нахождения матрично-векторных произведений при вычислении невязок матрича должна быть подвергнута декомпозиции на строчные блоки.

2. Нахождение скалярных произведений векторов (при построении ортонормального базиса подпространства Крылова).

Для наиболее эффективного вычисления скалярных проихведений и линейных комбинаций векоров базиса каждый из векторов [math]p_i[/math] должен быть разбит на блоки равного размера по числу процессоров.

Кроме того, на хранение векторов базиса приходится бОльшая часть затрат памяти.

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

Входные данные: Обязательные параметры:

[math]A[/math] - квадратная матрица (n x n).

[math]b[/math] - вектор правой части (n x 1).

Вариативно:

[math]x_0[/math] - начальное приближение (может не передаваться и быть всегда нулевым).

[math]\epsilon[/math] - погрешность (может не передаваться и быть заданным по умолчанию, обычно равна 1e-6).

Максимальное число итераций.

Выходные данные: [math]x^*[/math] - решение системы [math]Ax = b[/math]

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

1. Рассматриваемый алгоритм предназначен для решения невырожденных линейных систем большой размерности. При этом симметрия и положительная определенность матрицы не предполагается.

2. Для систем с разряженной матрицей алгоритм требует гораздо меньшей вычислительной мощности.

3. Особенностью метода является то, что с целью экономии памяти выбирается некоторая размерность Крылова <math>m<\math>, и после того как <math>m<\math> базисных векторов построены, решение уточняется.

4. Одной из важных практических особенностей метода является то, что на первых итерациях метод сходится гораздо быстрее, чем асимптотически.

Для ускорения сходимости метода минимальных невязок целесообразно время о времени использовать один шаг двух-шагового метода минимальных невязок.