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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[math]\displaystyle{ 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]\displaystyle{ H_i }[/math] - верхняя хессенбергова матрица порядка [math]\displaystyle{ i }[/math].

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

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

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

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

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

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

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

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

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

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

Пример

[math]\displaystyle{ \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]\displaystyle{ x_0 = 0 }[/math], то [math]\displaystyle{ r_0 = b }[/math]. Подпространства Крылова получаются такие:

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

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

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

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

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

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

[math]\displaystyle{ \|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]\displaystyle{ \lambda_{\mathrm{min}}(M) }[/math] и [math]\displaystyle{ \lambda_{\mathrm{max}}(M) }[/math] обозначают наибольшее и наименьшее собственное значение матрицы [math]\displaystyle{ M }[/math], соответственно.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Вариативно:

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

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

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

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

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

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