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

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

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

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

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

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

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

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

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

$r_i \perp A\mathcal{K}_i$.

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

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

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

$q_{i+1} \in \mathcal{K}_{i+1}, p_{i+1} = Aq_{i+1} \perp A\mathcal{K}_i$.

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

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

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

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

$AQ_i = Q_{i+1}\hat{H_i}, \hat{H_i}= \begin{bmatrix} H_i \\ 0 ... 0 & h_{i+1 i} \end{bmatrix},$

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

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

$\hat{H_i} = U_iR_i, R_i \in \mathbb{C}^{i*i}$.

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

$R_iy_i = z_i \equiv ||r_0||_2U^*_ie_1, e_1 = \begin{bmatrix} 1 0 ... 0 \end{bmatrix}^T,$.

Матрица $R_i$ невырожденная, следовательно,

$x_i = x_0 + Q_iy_i = x_0 + Q_iR_i^{-1}z_i$.

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

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

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

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

Пример

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

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

$\mathcal{K}_i = span \{ e_n, e_{n-1}, ..., e_{n-i+1}\}, i \le i \le n$

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

$x_0 = x_1 = ... = x_{n-1} = 0, x_n = e_1$.

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

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

Если матрица $A$ положительно определена, то

$\|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\|,$,

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

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

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

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

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

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

$\|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\|, \,$

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

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

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

Главные затраты $i$-ой итерации связаны с ортогонализацией.