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

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

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

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

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

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

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

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

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

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

Из рассуждений выше понятно что оптимальным является направление [math]\displaystyle{ p_i = - \bigtriangledown \phi (x_{i}) = b - Ax_i = r_i }[/math]. Величина [math]\displaystyle{ \alpha_i }[/math] выбирается из соображений [math]\displaystyle{ \alpha_i = \underset{\alpha}{\operatorname{argmin}} \phi (x_i + \alpha p_i) }[/math]. Аналитическую формулу [math]\displaystyle{ \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} }[/math] можно получить из [math]\displaystyle{ \frac{d}{d\alpha} \phi (x_i + \alpha p_i) = 0 }[/math].

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

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

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

Таким образом получим [math]\displaystyle{ x_{n} = x^* }[/math].
Кроме того, стоит обратить внимание на следующий факт: [math]\displaystyle{ 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 }[/math], тем самым формула для [math]\displaystyle{ \alpha_k }[/math] в этом методе совпадает с формулой в методе градиентного спуска.

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

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

1. Выбираем начальное приближение [math]\displaystyle{ x_0 }[/math]. Вычисляем невязку и вектор направления [math]\displaystyle{ r_0 = b - Ax_0, p_0 = r_0 }[/math].
2. Для [math]\displaystyle{ k = 0, 1, 2 ... }[/math] и до тех пор пока итерационный процесс не сойдется:
   1. Вычисляем [math]\displaystyle{ \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} }[/math]
   2. Новое приближение [math]\displaystyle{ x_{k + 1} = x_{k} + \alpha_k p_k }[/math]
   3. Новая невязка [math]\displaystyle{ r_{k + 1} = r_{k} - \alpha_k A p_k }[/math]
   4. Вычисляем [math]\displaystyle{ \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} }[/math]
   5. Новый вектор направления [math]\displaystyle{ p_k = r_{k+1} + \beta_k p_k }[/math]

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

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

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

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

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

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

Литература

• Методы численного анализа Евгений Евгеньевич Тыртышников
• 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