Участник:Сорокин Александр/Метод сопряженных градиентов (Решение СЛАУ)

Метод сопряженных градиентов для решения СЛАУ. Сорокин Александр 403

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

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

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

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

Пусть необходимо найти решение системы уравнений $Ax = b$, где $A^* = A \gt 0$.
Рассмотрим функционал $\phi (x) = \frac{1}{2}x^T A x - x^T b$.
Если $x^*$ это решение задачи минимизации данного функционала, то в этой точке градиент $\bigtriangledown \phi (x^*) = Ax^* - b$ должен быть равен нулю. Таким образом, минимизируя функционал $\phi (x)$ мы получим решение исходной системы.

1.2.1 Метод градиентного спуска

Как известно, градиент $\bigtriangledown \phi (x)$ является направлением наибольшего роста функции.
Метод градиентного спуска основан на стратегии движения в строну, противоположную возрастанию функционала. Оптимальным направлением в этом случае будет антиградиент $-\bigtriangledown \phi (x)$ и двигаться по нему нужно будет до тех пор, пока функционал убывает.
Таким образом можно построить следующий итерационный метод:

1. Выберем произвольное начальное приближение $x_0$.
2. $x_{i+1} = x_{i} + \alpha_i p_i$, где $p_i$ — направление движения, а $\alpha_i$ — величина шага.
   

Из рассуждений выше понятно что оптимальным является направление $p_i = - \bigtriangledown \phi (x_{i}) = b - Ax_i = r_i$. Величина $\alpha_i$ выбирается из соображений $\alpha_i = \underset{\alpha}{\operatorname{argmin}} \phi (x_i + \alpha p_i)$. Аналитическую формулу $\alpha_i = \frac{\bigtriangledown\phi (x_i)^T \bigtriangledown\phi (x_i)}{\bigtriangledown\phi (x_i)^T A \bigtriangledown\phi (x_i)} = \frac{r_i ^T r_i}{r_i ^T A r_i} = \frac{p_i ^T r_i}{p_i ^T A p_i}$ можно получить из $\frac{d}{d\alpha} \phi (x_i + \alpha p_i) = 0$.

1.2.2 Метод сопряженных направлений

Метод градиентного спуска обычно сходится очень долго. Можно построить алгоритм который сходится не больше чем за n шагов.
Предположим что у нас есть n линейно-независимых векторов $p_1, ... p_n$ таких что $(p_i, p_j)_A = (Ap_i, p_j) = 0, i \neq j$.
Так как имеется n таких векторов, то они образуют базис пространства и любой вектор можно выразить через них, в том числе $x^* - x_0 = \sum_{i = 1}^n \alpha_i p_i$.
Домножим это равенство А-скалярно на $p_k$ для всех k. С левой стороны получим $(x^* - x_0, p_k)_A = (b - Ax_0, p_k) = p_k ^T r_0$. С правой: $\sum_{i = 1}^n \alpha_i (p_i, p_k)_A = \alpha_k (p_k, p_k)_A = \alpha_k p_k ^T A p_k$. Тем самым $\alpha_k = \frac{p_k^T r_0}{p_k^T A p_k}$.
Итак мы получили новый алгоритм решения:

1. Выбираем начальное приближение $x_0$ и А-ортогональные направления $p_1, ..., p_n$.
2. Для всех $k = 1, ... n$:
   1. Выбираем $\alpha_k = \frac{p_k^T r_0}{p_k^T A p_k}$.
   2. Обновляем $x_{k+1} = x_{k} + \alpha_k p_k$.

Таким образом получим $x_{n} = x^*$.
Кроме того, стоит обратить внимание на следующий факт: $p_k^T r_0 = p_k^T(b - Ax_0) = p_k^T(b - A[x_0 + \alpha_1 p_1 + ... + \alpha_{k-1} p_{k-1}]) = p_k^T r_k$, тем самым формула для $\alpha_k$ в этом методе совпадает с формулой в методе градиентного спуска.

1.2.3 Метод сопряженных градиентов

Метод сопряженных градиентов это частный случай метода сопряженных направлений.
В этом методе мы будем искать новое приближение для решения $x$ и новый вектор направления $p$. Новый вектор направления будем получать используя алгоритм Грамма-Шмидта, где в качестве начальной линейно-независимой системы будем использовать невязки $r_i$. Метод будет выглядеть так:

1. Выбираем начальное приближение $x_0$. Вычисляем невязку и вектор направления $r_0 = b - Ax_0, p_0 = r_0$.
2. Для $k = 0, 1, 2 ...$ и до тех пор пока итерационный процесс не сойдется:
   1. Вычисляем $\alpha_k = \frac{p_k^T r_k}{p_k^T A p_k} = \frac{r_k^T r_k}{p_k^T A p_k}$
   2. Новое приближение $x_{k + 1} = x_{k} + \alpha_k p_k$
   3. Новая невязка $r_{k + 1} = r_{k} - \alpha_k A p_k$
   4. Вычисляем $\beta_k = -\frac{p_k^T A r_{k+1}}{p_k^T A p_k} = \frac{r_{k + 1}^T r_{k+1}}{r_k^T r_k}$
   5. Новый вектор направления $p_k = r_{k+1} + \beta_k p_k$

Последние два пункта алгоритма описывают процесс ортогонализации. Соотношения получились «короткими» в силу того, что $p_j^T Ar_{k+1} = 0, j \lt k$.

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

Вычислительное ядро алгоритма полностью совпадает с описанным выше итерационным процессом.

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

На каждой итерации алгоритма вычисляются:

1. 3 cкалярных произведения векторов
2. 3 обновление вектора
3. 1 умножение матрицы на вектор

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

Вычислить $r_0 = Ax_0$ для некоторого начального $x_0$.
$p_0 = r_0$
Для $k = 0, 1, 2 ...$
$\alpha_k = (r_k, r_k) / (p_k, Ap_k)$
$x_{k + 1} = x_k + \alpha_k p_k$
$r_{k + 1} = r_k - \alpha_k A p_k$
$\beta_k = (r_{k+1}, r_{k+1}) / (r_k, r_k)$
$p_{k+1} = r_{k+1} + \beta_k p_k$
Проверить условие сходимости; продолжить если необходимо
Конец

3 Литература

• Методы численного анализа Евгений Евгеньевич Тыртышников
• An Introduction to the Conjugate Gradient Method Without the Agonizing Pain by Jonathan Richard Shewchuk
• Iterative methods for sparse linear systems by Yousef Saad
• Notes on Some Methods for Solving Linear Systems by Dianne P. O'Leary