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

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{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$.