Методы подпространства Крылова
Метод сопряженных градиентов
Еще один класс итерационных методов может быть построен при рассмотрении задачи минимизации вектора невязки. Действительно, пусть $\bm{x^*}$ является решением СЛАУ $\bm{A} \bm{x} = \bm{b}$, где $\bm{A} \in \mathbb{R}^{n \times n}$. Тогда верным является утверждение:
\begin{equation}
\label{eq:min_res_raw}
\bm{x^*} = \argmin_{\bm{x}} \bm{r}^T \bm{r} = \argmin_{\bm{x}} (\bm{b} - \bm{A}\bm{x})^T (\bm{b} - \bm{A}\bm{x}).
\end{equation}
Рассмотрим случай положительно определенной и, следовательно, симметричной матрицы $\bm{A}$. Легко убедиться, что для симметричной матрицы $\bm{A}$ верным является утверждение:
\begin{equation}
\langle \bm{A} \bm{x}, \bm{y} \rangle = \langle \bm{x}, \bm{A} \bm{y} \rangle,
\end{equation}
где $\langle \cdot, \cdot \rangle$ – скалярное произведение векторов, $\bm{x}, \bm{y} \in \mathbb{R}^n$ – произвольные вектора. Тогда задача минимизации \eqref{eq:min_res_raw} может быть записана следующим образом:
\begin{equation}
\begin{split}
%\label{eq:min_res_raw}
\argmin_{\bm{x}} \left[\langle \bm{b}, \bm{b} \rangle - 2 \langle \bm{A} \bm{x}, \bm{b} \rangle + \langle \bm{A} \bm{x}, \bm{A} \bm{x} \rangle \right] &= \argmin_{\bm{x}} \left[\langle \bm{A} \bm{x}, \bm{A} \bm{x} \rangle - 2 \langle \bm{A} \bm{x}, \bm{b} \rangle \right] \\
&= \argmin_{\bm{x}} \left[\langle \bm{x}, \bm{A} \bm{x} \rangle - 2 \langle \bm{x}, \bm{b} \rangle \right].
\end{split}
\end{equation}
Обозначив целевую функцию как $g(\bm{x}) = \langle \bm{x}, \bm{A} \bm{x} \rangle - 2 \langle \bm{x}, \bm{b} \rangle$, мы имеем конечную форму задачи минимизации:
\begin{equation}
\label{eq:min_res}
\bm{x^*} = \argmin_{\bm{x}} g(\bm{x}).
\end{equation}
Для нахождения минимума, начав с некоторого приближения $\bm{x}$, необходимо определить направление поиска $\bm{v}$ и шаг $t$. Рассмотрим подобный поиск в контексте функции $g(\bm{x})$:
\begin{equation}
\begin{split}
g(\bm{x} + t \bm{v}) &= \langle \bm{x} + t \bm{v}, \bm{A} (\bm{x} + t \bm{v}) \rangle - 2 \langle \bm{x} + t \bm{v}, \bm{b} \rangle \\
&= \langle \bm{x}, \bm{A} \bm{x} \rangle + t \langle \bm{x}, \bm{A} \bm{v} \rangle + t \langle \bm{v}, \bm{A} \bm{x} \rangle + t^2 \langle \bm{v}, \bm{A} \bm{v} \rangle - 2 \langle \bm{x}, \bm{b} \rangle - 2 t \langle \bm{v}, \bm{b} \rangle \\
&= \langle \bm{x}, \bm{A} \bm{x} \rangle - 2 \langle \bm{x}, \bm{b} \rangle + 2 t \langle \bm{v}, \bm{A} \bm{x} \rangle + t^2 \langle \bm{v}, \bm{A} \bm{v} \rangle - 2 t \langle \bm{v}, \bm{b} \rangle \\
&= g(\bm{x}) + t^2 \langle \bm{v}, \bm{A} \bm{v} \rangle - 2 t \langle \bm{v}, \bm{b} - \bm{A} \bm{x} \rangle.
\end{split}
\end{equation}
В первую очередь найдем оптимальный шаг $t^{(opt)}$, минимизирующий функцию $g(\bm{x} + t \bm{v})$:
\begin{equation}
\begin{split}
&\frac{\partial g}{\partial t} = 2 t \langle \bm{v}, \bm{A} \bm{v} \rangle - 2 \langle \bm{v}, \bm{b} - \bm{A} \bm{x} \rangle = 0 \\
\Longrightarrow \quad &t^{(opt)} = \frac{\langle \bm{v}, \bm{b} - \bm{A} \bm{x} \rangle}{\langle \bm{v}, \bm{A} \bm{v} \rangle} \\
\Longrightarrow \quad &g(\bm{x} + t^{(opt)} \bm{v}) = g(\bm{x}) - \frac{\langle \bm{v}, \bm{b} - \bm{A} \bm{x} \rangle^2}{\langle \bm{v}, \bm{A} \bm{v} \rangle}.
\end{split}
\end{equation}
Из последнего выражения очевидно, что решение $\bm{x^*}$ минимизирует $g(\bm{x})$ для любого направления $\bm{v} \neq \bm{0}$. Теперь, задавшись начальным приближением $\bm{x}^{(0)}$ и начальным направлением поиска $\bm{v}^{(1)}$, мы можем построить следующий итерационный алгоритм:
\begin{align}
\label{eq:minim_iter_t_k}
t_{k} &= \frac{\langle \bm{v}^{(k)}, \bm{b} - \bm{A} \bm{x}^{(k-1)} \rangle}{\langle \bm{v}^{(k)}, \bm{A} \bm{v}^{(k)} \rangle}, \\
\label{eq:minim_iter_x_k}
\bm{x}^{(k)} &= \bm{x}^{(k-1)} + t_k \bm{v}^{(k)}.
\end{align}
Следующим шагом является генерация таких направлений поиска $\bm{v}^{(k)}$, что метод сходится достаточно быстро. Очевидным выбором является направление наискорейшего спуска:
\begin{equation}
- \frac{\partial g}{\partial \bm{x}} = - 2 (\bm{A} \bm{x} - \bm{b}) = 2\bm{r},
\end{equation}
то есть направлением наискорейшего спуска является вектор невязки. Тогда направление поиска определяется как:
\begin{equation}
\bm{v}^{(k+1)} = \bm{r}^{(k)},
\end{equation}
что задает
метод градиентного спуска для решения СЛАУ. Так как этот метод известен относительно медленной сходимостью, мы воспользуется альтернативным подходом, который базируется на построении таких $\bm{v}^{(i)}$, что
\begin{equation}
\langle \bm{v}^{(i)}, \bm{A} \bm{v}^{(j)} \rangle = 0, i \neq j.
\end{equation}
Такой вид ортогональности векторов мы будем называть $A$-ортогональностью, а полученную систему векторов $\{\bm{v}^{(1)}, \dots, \bm{v}^{(n)}\}$ $A$-ортогональной. Эта система задает базис пространства, так как является линейно независимой. Это легко проверить с помощью доказательства от обратного. Действительно, пусть можно выразить вектор $\bm{v}^{(i)}$ через линейную комбинацию остальных векторов:
\begin{equation}
\begin{split}
&\bm{v}^{(i)} = \sum_{\substack{j = 1 \\ j \neq i}} \alpha_j \bm{v}^{(j)} \\
\Longrightarrow \quad &\langle \bm{v}^{(i)}, \bm{A} \bm{v}^{(i)} \rangle = \sum_{\substack{j = 1 \\ j \neq i}} \alpha_j \langle \bm{v}^{(j)}, \bm{A} \bm{v}^{(i)} \rangle \\
\Longrightarrow \quad &\langle \bm{v}^{(i)}, \bm{A} \bm{v}^{(i)} \rangle = 0,
\end{split}
\end{equation}
что нарушает условие $A$-ортогональности, из чего следует, что система векторов действительно является линейно независимой. Результирующий итерационный метод, построенный на основе уравнений \eqref{eq:minim_iter_t_k} и \eqref{eq:minim_iter_x_k}, называется
методом сопряженных направлений. Следующая теорема доказывает, что такой метод дает точное решение в контексте арифметики с бесконечной точностью за $n$ итераций.
Пусть система векторов $\{\bm{v}^{(1)}, \dots, \bm{v}^{(n)}\}$ является $A$-ортогональной, а матрица $\bm{A}$ положительно определенной. Тогда $\bm{A} \bm{x}^{(n)} = \bm{b}$, где $\bm{x}^{(n)}$ определяется с помощью итерационного алгоритма ниже:
\begin{align}
t_{k} &= \frac{\langle \bm{v}^{(k)}, \bm{b} - \bm{A} \bm{x}^{(k-1)} \rangle}{\langle \bm{v}^{(k)}, \bm{A} \bm{v}^{(k)} \rangle}, \\
\bm{x}^{(k)} &= \bm{x}^{(k-1)} + t_k \bm{v}^{(k)}.
\end{align}
Раскроем итерации в выражении $\bm{A} \bm{x}^{(n)}$ вплоть до начального приближения $\bm{x}^{(0)}$:
\begin{equation}
\begin{split}
\bm{A} \bm{x}^{(n)} &= \bm{A} \left(\bm{x}^{(n-1)} + t_n \bm{v}^{(n)}\right), \\
&= \bm{A} \left(\bm{x}^{(n-2)} + t_{n-1} \bm{v}^{(n-1)} + t_n \bm{v}^{(n)}\right), \\
&= \dots \\
&= \bm{A} \left(\bm{x}^{(0)} + t_{1} \bm{v}^{(1)} + \dots + t_n \bm{v}^{(n)}\right),
\end{split}
\end{equation}
из чего следует:
\begin{equation}
\begin{split}
\bm{A} \bm{x}^{(n)} - \bm{b} = \bm{A} \bm{x}^{(0)} - \bm{b} + \bm{A} \left(t_{1} \bm{v}^{(1)} + \dots + t_n \bm{v}^{(n)}\right).
\end{split}
\end{equation}
Посчитаем скалярное произведение обоих сторон уравнения, домножив его на $\bm{v}^{(k)}$:
\begin{equation}
\label{eq:th_n_steps_1}
\begin{split}
\langle \bm{A} \bm{x}^{(n)} - \bm{b}, \bm{v}^{(k)} \rangle = \langle \bm{A} \bm{x}^{(0)} - \bm{b}, \bm{v}^{(k)} \rangle + t_k \langle \bm{v}^{(k)}, \bm{A} \bm{v}^{(k)} \rangle.
\end{split}
\end{equation}
Из определения $t_k$ можно получить выражение:
\begin{equation}
\begin{split}
t_{k} \langle \bm{v}^{(k)}, \bm{A} \bm{v}^{(k)} \rangle &= \langle \bm{v}^{(k)}, \bm{b} - \bm{A} \bm{x}^{(k-1)} \rangle \\
&= \langle \bm{v}^{(k)}, \bm{b} - \bm{A} \bm{x}^{(0)} + \bm{A} \bm{x}^{(0)} - \bm{A} \bm{x}^{(1)} + \bm{A} \bm{x}^{(1)} - \dots \\
&\quad - \bm{A} \bm{x}^{(k-2)} + \bm{A} \bm{x}^{(k-2)} - \bm{A} \bm{x}^{(k-1)} \rangle \\
&= \langle \bm{v}^{(k)}, \bm{b} - \bm{A} \bm{x}^{(0)}\rangle + \langle \bm{v}^{(k)}, \bm{A} \bm{x}^{(0)} - \bm{A} \bm{x}^{(1)} \rangle + \dots \\
&\quad + \langle \bm{v}^{(k)}, \bm{A} \bm{x}^{(k-2)} - \bm{A} \bm{x}^{(k-1)} \rangle.
\end{split}
\end{equation}
Заметим, что из $\bm{x}^{(i)} = \bm{x}^{(i-1)} + t_i \bm{v}^{(i)}$ следует $\bm{A} \bm{x}^{(i-1)} - \bm{A} \bm{x}^{(i)} = -t_i \bm{A} \bm{v}^{(i)}$. Тогда подстановка в выражение выше дает:
\begin{equation}
t_{k} \langle \bm{v}^{(k)}, \bm{A} \bm{v}^{(k)} \rangle = \langle \bm{v}^{(k)}, \bm{b} - \bm{A} \bm{x}^{(0)}\rangle,
\end{equation}
из чего после подстановки в \eqref{eq:th_n_steps_1} мы получаем:
\begin{equation}
\begin{split}
&\langle \bm{A} \bm{x}^{(n)} - \bm{b}, \bm{v}^{(k)} \rangle = \langle \bm{A} \bm{x}^{(0)} - \bm{b}, \bm{v}^{(k)} \rangle + \langle \bm{v}^{(k)}, \bm{b} - \bm{A} \bm{x}^{(0)}\rangle \\
\Longrightarrow \quad &\langle \bm{A} \bm{x}^{(n)} - \bm{b}, \bm{v}^{(k)} \rangle = 0 \quad \forall k = 1, \dots, n.
\end{split}
\end{equation}
Последнее выражение говорит о том, что вектор невязки на $n$-й итерации ортогонален всем векторам $\{\bm{v}^{(1)}, \dots, \bm{v}^{(n)}\}$, которые, будучи линейно независимыми, задают базис пространства. Это возможно тогда и только тогда, когда $\bm{A} \bm{x}^{(n)} - \bm{b} = 0$, что означает, что $\bm{x}^{(n)}$ является решением исходной СЛАУ.
Также можно доказать, что получаемые векторы невязки $\bm{r}^{(k)}$ будут ортогональны векторам $\bm{v}^{(i)}$, $i = 1, \dots, k$.
Метод сопряженных градиентов выбирает такие $\{\bm{v}^{(i)}\}$, что система, состоящая из векторов невязки $\{\bm{r}^{(i)}\}$, является ортогональной, т.е. $\langle \bm{r}^{(i)}, \bm{r}^{(j)} \rangle$ для $i \neq j$. Соответствующие формулы для $\{\bm{v}^{(i)}\}$ можно построить, адаптировав процесс Грама-Шмидта для $A$-ортогональности. Действительно, пусть начальное направление является направлением наискорейшего спуска, т.е. $\bm{v}^{(1)} = \bm{r}^{(0)}$. Предположим, что мы уже рассчитали приближения $\{\bm{x}^{(1)}, \dots, \bm{x}^{(k-1)}\}$ и направления $\{\bm{v}^{(1)}, \dots, \bm{v}^{(k-1)}\}$. Применим процесс Грама-Шмидт для нахождения $\bm{v}^{(k)}$:
\begin{equation}
\bm{v}^{(k)} = \bm{r}^{k-1} + s_{k-1} \bm{v}^{(k-1)},
\end{equation}
где коэффициент $s_{k-1}$ мы находим исходя из условия $A$-ортогональности $\langle \bm{v}^{(k)}, \bm{A} \bm{v}^{(k-1)} \rangle = 0$:
\begin{equation}
\begin{split}
&0 = \langle \bm{r}^{k-1}, \bm{A} \bm{v}^{(k-1)} \rangle + s_{k-1} \langle \bm{v}^{(k-1)}, \bm{A} \bm{v}^{(k-1)} \rangle \\
\Longrightarrow \quad &s_{k-1} = -\frac{\langle \bm{r}^{(k-1)}, \bm{A} \bm{v}^{(k-1)} \rangle}{\langle \bm{v}^{(k-1)}, \bm{A} \bm{v}^{(k-1)} \rangle}.
\end{split}
\end{equation}
Имея выражение для $\bm{v}^{(k)}$, переформулируем $t_k$:
\begin{equation}
\begin{split}
t_{k} &= \frac{\langle \bm{v}^{(k)}, \bm{r}^{(k-1)} \rangle}{\langle \bm{v}^{(k)}, \bm{A} \bm{v}^{(k)} \rangle} \\
&= \frac{\langle \bm{r}^{(k-1)}, \bm{r}^{(k-1)} \rangle}{\langle \bm{v}^{(k)}, \bm{A} \bm{v}^{(k)} \rangle} + s_{k-1} \frac{\langle \bm{v}^{(k-1)}, \bm{r}^{(k-1)} \rangle}{\langle \bm{v}^{(k)}, \bm{A} \bm{v}^{(k)} \rangle} \\
&= \frac{\langle \bm{r}^{(k-1)}, \bm{r}^{(k-1)} \rangle}{\langle \bm{v}^{(k)}, \bm{A} \bm{v}^{(k)} \rangle},
\end{split}
\label{eq:t_k_final}
\end{equation}
где мы использовали свойство $\langle \bm{v}^{(k-1)}, \bm{r}^{(k-1)} \rangle = 0$. Мы также можем записать рекурсивное выражение для $\bm{r}^{(k)}$:
\begin{equation}
\begin{split}
&\bm{x}^{(k)} = \bm{x}^{(k-1)} + t_k \bm{v}^{(k)} \\
\Longrightarrow \quad &\bm{A}\bm{x}^{(k)} - \bm{b} = \bm{A}\bm{x}^{(k-1)} - \bm{b} + t_k \bm{A} \bm{v}^{(k)} \\
\Longrightarrow \quad &\bm{r}^{(k)} = \bm{r}^{(k-1)} - t_k \bm{A} \bm{v}^{(k)}.
\end{split}
\end{equation}
Требование ортогональности $\bm{r}^{(k)}$ и $\bm{r}^{(k-1)}$ при этом дает
\begin{equation}
\langle \bm{v}^{(k)}, \bm{A} \bm{v}^{(k)} \rangle = \frac{1}{t_k} \langle \bm{r}^{(k-1)} \bm{r}^{(k-1)} \rangle.
\end{equation}
Это позволяет упросить выражение для $s_k$, если мы заметим что:
\begin{equation}
\langle \bm{v}^{(k)}, \bm{A} \bm{r}^{(k)} \rangle = \langle \bm{A} \bm{v}^{(k)}, \bm{r}^{(k)} \rangle = \frac{1}{t_k}\langle \bm{r}^{(k-1)} - \bm{r}^{(k)}, \bm{r}^{(k)} \rangle = -\frac{1}{t_k}\langle \bm{r}^{(k)}, \bm{r}^{(k)} \rangle,
\end{equation}
Тогда выражение для $s_k$ принимает вид:
\begin{equation}
s_k = \frac{\langle \bm{r}^{(k)}, \bm{r}^{(k)} \rangle}{\langle \bm{r}^{(k-1)}, \bm{r}^{(k-1)} \rangle}.
\end{equation}
Резюмируем метод сопряженных градиентов в виде следующих рекурсивных уравнений:
\begin{align}
t_k &= \frac{\langle \bm{r}^{(k-1)}, \bm{r}^{(k-1)} \rangle}{\langle \bm{v}^{(k)}, \bm{A} \bm{v}^{(k)} \rangle}, \\
\bm{x}^{(k)} &= \bm{x}^{(k-1)} + t_k \bm{v}^{(k)}, \\
\bm{r}^{(k)} &= \bm{r}^{(k-1)} - t_k \bm{A} \bm{v}^{(k)}, \\
s_k &= \frac{\langle \bm{r}^{(k)}, \bm{r}^{(k)} \rangle}{\langle \bm{r}^{(k-1)}, \bm{r}^{(k-1)} \rangle}, \\
\bm{v}^{(k+1)} &= \bm{r}^{(k)} + s_k \bm{v}^{(k)},
\end{align}
где $\bm{x}^{(0)}$, $\bm{r}^{(0)} = \bm{b} - \bm{A} \bm{x}^{(0)}$ и $\bm{v}^{(1)} = \bm{r}^{(0)}$ являются начальными условиями.
Важно отметить, что сопряженные направления напрямую связаны с подпространством Крылова.
Подпространством Крылова размерности $n$ называется линейное пространство $\mathcal{K}_m(\bm{z}, \bm{A})$, порожденное множеством векторов $\{\bm{z}, \bm{A}\bm{z}, \bm{A}^2\bm{z}, \dots, \bm{A}^{m-1}\bm{z}\}$, где $\bm{z} \in \mathbb{R}^n$ и $\bm{A} \in \mathbb{R}^{n \times n}$.
Рассмотрим выражения для сопряженных направлений:
\begin{align*}
\bm{v}^{(1)} &= \bm{r}^{(0)}, \\
\bm{v}^{(2)} &= -t_1 \bm{A}\bm{r}^{(0)} + (1 + s_1) \bm{r}^{(0)}, \\
\bm{v}^{(3)} &= \bm{r}^{(0)} - \left[t_1 + t_2 (1 + s_1)\right] \bm{A}\bm{r}^{(0)} + t_1 t_2 \bm{A}^2 \bm{r}^{(0)}, \\
&\dots
\end{align*}
Можно заметить, что сопряженные направления $\bm{v}^{k}$ принадлежат подпространству Крылова $\mathcal{K}_n(\bm{r}^{(0)}, \bm{A})$. Более того, несложно продемонстрировать, что они составляют ортонормальный базис подпространства Крылова $\mathcal{K}_n(\bm{r}^{(0)}, \bm{A})$.
Предобуславливание матриц
В случае, когда матрица $\bm{A}$ является плохо обусловленной (т.е. имеет большое число обусловленности $K(\bm{A})$), вычислительная погрешность, как мы выяснили, становится достаточно большой и значимой. Более того, плохо обусловленные матрицы обладают медленной сходимостью. Выходом из этой ситуации является модификация исходной матрицы $\bm{A}$ так, что полученная в результате матрица обладает значительно меньшим числом обусловленности. Матрица, используемая для модификации исходной матрицы, называется
матрицей предобуславливания. Для положительно определенных матрицы такая модификация должна иметь вид:
\begin{equation}
\tilde{\bm{A}} = \bm{C}^{-1} \bm{A} \left(\bm{C}^{-1}\right)^{T},
\end{equation}
где обозначение $\bm{C}^{-1}$ говорит о том, что
матрица предобуславливания $\bm{C}^{-1}$ должна быть невырожденной. При подобном предобуславливании полученная матрица $\tilde{\bm{A}}$ сохраняет свойство положительной определенности. Одним из вариантов матриц $\bm{C}^{-1}$ является следующее предобуславливание:
\begin{equation}
\tilde{\bm{A}} = \bm{D}^{-1/2} \bm{A} \left(\bm{D}^{-1/2}\right)^{T},
\end{equation}
где $\bm{D}^{-1/2}$ состоит из корней от обратных диагональных элементов матрицы $\bm{A}$. Метод сопряженных градиентов, сформулированный для модифицированной матрицы, называется
методом сопряженных градиентов с предобуславливанием.
Несложно продемонстрировать, что его формулировка будет иметь следующий вид:
Резюмируем метод сопряженных градиентов в виде следующих рекурсивных уравнений:
\begin{align}
t_k &= \frac{\langle \bm{w}^{(k-1)}, \bm{w}^{(k-1)} \rangle}{\langle \bm{v}^{(k)}, \bm{A} \bm{v}^{(k)} \rangle}, \\
\bm{x}^{(k)} &= \bm{x}^{(k-1)} + t_k \bm{v}^{(k)}, \\
\bm{r}^{(k)} &= \bm{r}^{(k-1)} - t_k \bm{A} \bm{v}^{(k)}, \\
\bm{w}^{(k)} &= \bm{C}^{-1} \bm{r}^{(k)}, \\
s_k &= \frac{\langle \bm{w}^{(k)}, \bm{w}^{(k)} \rangle}{\langle \bm{w}^{(k-1)}, \bm{w}^{(k-1)} \rangle}, \\
\bm{v}^{(k+1)} &= \bm{C}^{-T} \bm{w}^{(k)} + s_k \bm{v}^{(k)},
\end{align}
где $\bm{x}^{(0)}$, $\bm{r}^{(0)} = \bm{b} - \bm{A} \bm{x}^{(0)}$, $\bm{w}^{(0)} = \bm{C}^{-1} \bm{r}^{(0)}$ и $\bm{v}^{(1)} = \bm{C}^{-T} \bm{w}^{(0)}$ являются начальными условиями.