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

В числовой линейной алгебре метод сопряженных градиентов — это итерационный метод численного решения линейной системы.

${\boldsymbol {Ax}}={\boldsymbol {b}}$

где ${\boldsymbol {A}}$ является симметричным положительно определенным без вычисления ${\boldsymbol {A}}^{-1}$ явно. Метод сопряженных градиентов можно использовать с нескольких разных точек зрения, включая специализацию метода сопряженных направлений. для оптимизации и вариации итерации Арнольди / Ланцоша для собственных значений задач .

Цель этой статьи — задокументировать важные этапы этих выводов.

Сопряженное направление

Метод сопряженных градиентов можно рассматривать как частный случай метода сопряженных направлений, применяемого для минимизации квадратичной функции.

$f({\boldsymbol {x}})={\boldsymbol {x}}^{\mathrm {T} }{\boldsymbol {A}}{\boldsymbol {x}}-2{\boldsymbol {b}}^{\mathrm {T} }{\boldsymbol {x}}{\text{.}}$

что позволяет нам применить геометрическую интуицию.

Поиск линии

Геометрически квадратичную функцию можно эквивалентным образом представить, записав ее значение в каждой точке пространства. Точки равного значения составляют его контурные поверхности, представляющие собой концентрические эллипсоиды с уравнением $${\boldsymbol {x}}^{\mathrm {T} }{\boldsymbol {A}}{\boldsymbol {x}}-2{\boldsymbol {b}}^{\mathrm {T} }{\boldsymbol {x}}=C$$для изменения $C$. Как $C$ уменьшается, эллипсоиды становятся все меньше и меньше, пока при минимальном значении эллипсоид не сожмется до их общего центра.

Минимизация квадратичной функции тогда представляет собой проблему перемещения по плоскости в поисках общего центра всех этих эллипсоидов. Центр можно найти, вычислив ${\boldsymbol {A}}^{-1}$ явно, но именно этого мы и пытаемся избежать.

Самый простой метод — жадный поиск по строке , при котором мы начинаем с некоторой точки ${\boldsymbol {x}}_{0}$, выбери направление ${\boldsymbol {p}}_{0}$ как-нибудь, а потом свести к минимуму $f({\boldsymbol {x}}_{0}+{\boldsymbol {p}}_{0}\alpha _{0})$. Это имеет простое решение в замкнутой форме, которое не требует обращения матрицы: $$\alpha _{0}={\frac {{\boldsymbol {p}}_{0}^{\mathrm {T} }({\boldsymbol {b}}-{\boldsymbol {Ax}}_{0})}{{\boldsymbol {p}}_{0}^{\mathrm {T} }{\boldsymbol {A}}{\boldsymbol {p}}_{0}}}$$Геометрически мы начинаем в какой-то момент ${\boldsymbol {x}}_{0}$ на некотором эллипсоиде, затем выбираем направление и идем в этом направлении, пока не достигнем точки, в которой эллипсоид сворачивается в этом направлении. Это не обязательно минимум, но это прогресс к нему. Визуально он движется вдоль прямой и останавливается, как только мы достигаем точки, касательной к контуру эллипсоида.

Теперь мы можем повторить эту процедуру, начиная с новой точки. ${\boldsymbol {x}}_{1}={\boldsymbol {x}}_{0}+\alpha _{0}{\boldsymbol {p}}_{0}$, выберите новое направление ${\boldsymbol {p}}_{1}$, вычислить $\alpha _{1}$, и т. д.

Мы можем обобщить это в виде следующего алгоритма:

Начните с выбора первоначального предположения ${\boldsymbol {x}}_{0}$и вычислим начальную невязку ${\boldsymbol {r}}_{0}={\boldsymbol {b}}-{\boldsymbol {Ax}}_{0}$, затем повторите:

${\begin{aligned}\alpha _{i}&={\frac {{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {r}}_{i}}{{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ap}}_{i}}}{\text{,}}\\{\boldsymbol {x}}_{i+1}&={\boldsymbol {x}}_{i}+\alpha _{i}{\boldsymbol {p}}_{i}{\text{,}}\\{\boldsymbol {r}}_{i+1}&={\boldsymbol {r}}_{i}-\alpha _{i}{\boldsymbol {Ap}}_{i}\end{aligned}}$

где ${\boldsymbol {p}}_{0},{\boldsymbol {p}}_{1},{\boldsymbol {p}}_{2},\ldots$ должны быть выбраны. Обратите внимание, в частности, как остаток вычисляется итеративно, шаг за шагом, а не каждый раз заново: $${\boldsymbol {r}}_{i+1}={\boldsymbol {b}}-{\boldsymbol {Ax}}_{i+1}={\boldsymbol {b}}-{\boldsymbol {A}}({\boldsymbol {x}}_{i}+\alpha _{i}{\boldsymbol {p}}_{i})={\boldsymbol {r}}_{i}-\alpha _{i}{\boldsymbol {A}}{\boldsymbol {p}}_{i}$$Возможно, это правда, что $\alpha _{i}=0$ преждевременно, что приведет к численным проблемам. Однако для конкретного выбора ${\boldsymbol {p}}_{0},{\boldsymbol {p}}_{1},{\boldsymbol {p}}_{2},\ldots$, то до сходимости этого не произойдет, как мы докажем ниже.

Сопряженные направления

Если направления ${\boldsymbol {p}}_{0},{\boldsymbol {p}}_{1},{\boldsymbol {p}}_{2},\ldots$ выбраны неудачно, то прогресс будет медленным. В частности, метод градиентного спуска будет медленным. Это можно увидеть на диаграмме, где зеленая линия — результат выбора локального направления градиента. Он движется зигзагами к минимуму, но неоднократно превышает его. Напротив, если мы выберем направления как набор взаимно сопряженных направлений , то перерегулирования не будет, и мы получим глобальный минимум после $n$ шаги, где $n$ это количество измерений.

Два сопряженных диаметра эллипса . Каждое ребро ограничивающего параллелограмма параллельно . одному из диаметров

Понятие сопряженных направлений пришло из классической геометрии эллипса. Для эллипса центры двух полуосей взаимно сопряжены относительно эллипса тогда и только тогда, когда линии параллельны касательной, ограничивающей параллелограмм, как показано на рисунке. Концепция обобщается на n -мерные эллипсоиды, где n полуосей $t_{0}{\boldsymbol {p}}_{0},\dots ,t_{n-1}{\boldsymbol {p}}_{n-1}$ взаимно сопряжены относительно эллипсоида тогда и только тогда, когда каждая ось параллельна касательной, ограничивающей параллелепипед . Другими словами, для любого $i$, касательная плоскость к эллипсоиду в точке ${\boldsymbol {c}}+t_{i}{\boldsymbol {p}}_{i}$ представляет собой гиперплоскость, натянутую на векторы $\{{\boldsymbol {p}}_{j}:j\neq i\}$, где ${\boldsymbol {c}}$ является центром эллипсоида.

Обратите внимание, что нам нужно масштабировать каждый вектор направления. ${\boldsymbol {p}}_{i}$ по скаляру $t_{i}$, так что ${\boldsymbol {c}}+t_{i}{\boldsymbol {p}}_{i}$ падает точно на эллипсоид.

Дан эллипсоид с уравнением $${\boldsymbol {x}}^{\mathrm {T} }{\boldsymbol {A}}{\boldsymbol {x}}-2{\boldsymbol {b}}^{\mathrm {T} }{\boldsymbol {x}}=C$$для некоторой константы $C$, мы можем перевести его так, чтобы его центр находился в начале координат. Это меняет уравнение на $${\boldsymbol {x}}^{\mathrm {T} }{\boldsymbol {A}}{\boldsymbol {x}}=C'$$для какой-то другой константы $C'$. Тогда условие касания: $$(t_{i}{\boldsymbol {p}}_{i}+{\boldsymbol {p}}_{j}dt_{j})^{\mathrm {T} }{\boldsymbol {A}}(t_{i}{\boldsymbol {p}}_{i}+{\boldsymbol {p}}_{j}dt_{j})=C'+O(dt_{j}^{2}),\quad \forall i\neq j$$то есть, ${\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ap}}_{j}=0$ для любого $i\neq j$.

Метод сопряженных направлений неточен в том смысле, что не приводятся формулы для выбора направлений. ${\boldsymbol {p}}_{0},{\boldsymbol {p}}_{1},{\boldsymbol {p}}_{2},\ldots$. Конкретный выбор приводит к использованию различных методов, включая метод сопряженных градиентов и метод исключения Гаусса .

Процесс Грама – Шмидта

Мы можем свести в таблицу уравнения, которые нам нужно обнулить:

	0	1	2	3	$\cdots$
0		${\boldsymbol {p}}_{0}^{\mathrm {T} }{\boldsymbol {Ap}}_{1}$	${\boldsymbol {p}}_{0}^{\mathrm {T} }{\boldsymbol {Ap}}_{2}$	${\boldsymbol {p}}_{0}^{\mathrm {T} }{\boldsymbol {Ap}}_{3}$	$\cdots$
1			${\boldsymbol {p}}_{1}^{\mathrm {T} }{\boldsymbol {Ap}}_{2}$	${\boldsymbol {p}}_{1}^{\mathrm {T} }{\boldsymbol {Ap}}_{3}$	$\cdots$
2				${\boldsymbol {p}}_{2}^{\mathrm {T} }{\boldsymbol {Ap}}_{3}$	$\cdots$
$\vdots$					$\ddots$

Это напоминает проблему ортогонализации, которая требует ${\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {p}}_{j}=0$ для любого $i\neq j$, и ${\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {p}}_{j}=1$ для любого $i=j$. Таким образом, проблема поиска сопряженных осей менее ограничена, чем проблема ортогонализации, поэтому процесс Грама – Шмидта работает с дополнительными степенями свободы, которые мы позже можем использовать, чтобы выбрать те, которые упростят вычисления:

• Произвольно задано ${\boldsymbol {p}}_{0}$.
• Произвольно задано ${\boldsymbol {p}}_{10}$, затем измените его на ${\boldsymbol {p}}_{1}={\boldsymbol {p}}_{10}-{\frac {{\boldsymbol {p}}_{0}^{\mathrm {T} }{\boldsymbol {Ap}}_{10}}{{\boldsymbol {p}}_{0}^{\mathrm {T} }{\boldsymbol {Ap}}_{0}}}{\boldsymbol {p}}_{0}$.
• Произвольно задано ${\boldsymbol {p}}_{20}$, затем измените его на ${\boldsymbol {p}}_{2}={\boldsymbol {p}}_{20}-\sum _{i=0}^{1}{\frac {{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ap}}_{20}}{{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ap}}_{i}}}{\boldsymbol {p}}_{i}$.
• ...
• Произвольно задано ${\boldsymbol {p}}_{n-1,0}$, затем измените его на ${\boldsymbol {p}}_{n-1}={\boldsymbol {p}}_{n-1,0}-\sum _{i=0}^{n-2}{\frac {{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ap}}_{n-1,0}}{{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ap}}_{i}}}{\boldsymbol {p}}_{i}$.

Самый естественный выбор ${\boldsymbol {p}}_{k,0}$ это градиент. То есть, ${\boldsymbol {p}}_{k,0}=\nabla f({\boldsymbol {x}}_{k})$. Поскольку сопряженные направления можно масштабировать на ненулевое значение, мы масштабируем его на $-1/2$ для чистоты обозначений, получив $${\boldsymbol {p}}_{k,0}=\mathbf {r} _{k}=\mathbf {b} -\mathbf {Ax} _{k}$$Таким образом, мы имеем ${\boldsymbol {p}}_{k}={\boldsymbol {r}}_{k}-\sum _{i=0}^{k-1}{\frac {{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ar}}_{k}}{{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ap}}_{i}}}{\boldsymbol {p}}_{i}$. Подключив его, мы имеем алгоритм сопряженного градиента: $${\begin{aligned}&\mathbf {r} _{0}:=\mathbf {b} -\mathbf {Ax} _{0}\\&\mathbf {p} _{0}:=\mathbf {r} _{0}\\&k:=0\\&{\text{do while }}k<n\\&\qquad \alpha _{k}:={\frac {\mathbf {p} _{k}^{\mathsf {T}}\mathbf {r} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {Ap} _{k}}}\\&\qquad \mathbf {x} _{k+1}:=\mathbf {x} _{k}+\alpha _{k}\mathbf {p} _{k}\\&\qquad {\text{if }}|\alpha _{k}|{\text{ is sufficiently small, then exit loop}}\\&\qquad \mathbf {r} _{k+1}:=\mathbf {r} _{k}-\alpha _{k}\mathbf {Ap} _{k}\\&\qquad \mathbf {p} _{k+1}:={\boldsymbol {r}}_{k+1}-\sum _{i=0}^{k}{\frac {{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ar}}_{k+1}}{{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ap}}_{i}}}{\boldsymbol {p}}_{i}\\&\qquad k:=k+1\\&{\text{return }}\mathbf {x} _{k+1}{\text{ as the result}}\end{aligned}}$$Предложение. Если в какой-то момент, $\alpha _{k}=0$, то алгоритм сошёлся, т.е. $\nabla f(\mathrm {x} _{k+1})=0$.

Доказательство. По построению это будет означать, что $\mathbf {x} _{k+1}=\mathbf {x} _{k}$, то есть шаг по сопряженному градиенту возвращает нас точно туда, где мы были. Это возможно только в том случае, если локальный градиент уже равен нулю.

Упрощение

Этот алгоритм можно значительно упростить с помощью некоторых лемм, в результате чего получится алгоритм сопряженного градиента.

Лемма 1. $\mathbf {p} _{i}^{T}\mathbf {r} _{j}=0,\;\forall i<j$ и $\mathbf {r} _{i}^{T}\mathbf {r} _{j}=0,\;\forall i<j$.

Доказательство. По геометрической конструкции касательная плоскость к эллипсоиду в точке $\mathbf {x} _{j}$ содержит каждый из предыдущих сопряженных векторов направления $\mathbf {p} _{0},\mathbf {p} _{1},\dots ,\mathbf {p} _{j-1}$. Дальше, $\mathbf {r} _{j}$ перпендикулярен касательной, поэтому $\mathbf {p} _{i}^{T}\mathbf {r} _{j}=0,\;\forall i<j$. Второе уравнение верно, поскольку по построению $\mathbf {r} _{0},\mathbf {r} _{1},\dots ,\mathbf {r} _{j-1}$ является линейным преобразованием $\mathbf {p} _{0},\mathbf {p} _{1},\dots ,\mathbf {p} _{j-1}$.

Лемма 2. $\mathbf {p} _{k}^{T}\mathbf {r} _{k}=\mathbf {r} _{k}^{T}\mathbf {r} _{k}$.

Доказательство. По конструкции, $\mathbf {p} _{k}:={\boldsymbol {r}}_{k}-\sum _{i=0}^{k-1}{\frac {{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ar}}_{k-1}}{{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ap}}_{i}}}{\boldsymbol {p}}_{i}$, теперь применим лемму 1.

Лемма 3. ${\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ar}}_{k+1}={\begin{cases}0,\;i<k\\-{\boldsymbol {r}}_{k+1}^{T}{\boldsymbol {r}}_{k+1}/\alpha _{k},\;i=k\end{cases}}$.

Доказательство. По конструкции мы имеем $\mathbf {r} _{i+1}=\mathbf {r} _{i}-\alpha _{k}\mathbf {Ap} _{i}$, таким образом $${\boldsymbol {r}}_{k+1}^{T}{\boldsymbol {A}}{\boldsymbol {p}}_{i}={\boldsymbol {r}}_{k+1}^{T}{\frac {{\boldsymbol {r}}_{i}-{\boldsymbol {r}}_{i+1}}{\alpha _{i}}}$$Теперь применим лемму 1.

Подставив леммы 1-3, получим $\alpha _{k}={\frac {\mathbf {r} _{k}^{\top }\mathbf {r} _{k}}{\mathbf {p} _{k}^{\top }\mathbf {A} \mathbf {p} _{k}}}$ и $\mathbf {p} _{k+1}:={\boldsymbol {r}}_{k+1}+{\frac {\mathbf {r} _{k+1}^{\top }\mathbf {r} _{k+1}}{\mathbf {r} _{k}^{\top }\mathbf {r} _{k}}}\mathbf {p} _{k}$, который является правильным алгоритмом сопряженного градиента.

Итерация Арнольди/Ланцоша

Дополнительная информация: итерация Арнольди и итерация Ланцоша.

Метод сопряженных градиентов также можно рассматривать как вариант итерации Арнольди/Ланцоша, применяемый для решения линейных систем.

Общий метод Арнольда

В итерации Арнольди все начинается с вектора ${\boldsymbol {r}}_{0}$ и постепенно строит ортонормированный базис $\{{\boldsymbol {v}}_{1},{\boldsymbol {v}}_{2},{\boldsymbol {v}}_{3},\ldots \}$ of the Krylov subspace

${\mathcal {K}}({\boldsymbol {A}},{\boldsymbol {r}}_{0})=\mathrm {span} \{{\boldsymbol {r}}_{0},{\boldsymbol {Ar}}_{0},{\boldsymbol {A}}^{2}{\boldsymbol {r}}_{0},\ldots \}$

определяя ${\boldsymbol {v}}_{i}={\boldsymbol {w}}_{i}/\lVert {\boldsymbol {w}}_{i}\rVert _{2}$ где

${\boldsymbol {v}}_{i}={\begin{cases}{\boldsymbol {r}}_{0}&{\text{if }}i=1{\text{,}}\\{\boldsymbol {Av}}_{i-1}-\sum _{j=1}^{i-1}({\boldsymbol {v}}_{j}^{\mathrm {T} }{\boldsymbol {Av}}_{i-1}){\boldsymbol {v}}_{j}&{\text{if }}i>1{\text{.}}\end{cases}}$

Другими словами, для $i>1$, ${\boldsymbol {v}}_{i}$ находится методом ортогонализации Грама-Шмидта ${\boldsymbol {Av}}_{i-1}$ против $\{{\boldsymbol {v}}_{1},{\boldsymbol {v}}_{2},\ldots ,{\boldsymbol {v}}_{i-1}\}$ с последующей нормализацией.

В матричной форме итерация описывается уравнением

${\boldsymbol {AV}}_{i}={\boldsymbol {V}}_{i+1}{\boldsymbol {\tilde {H}}}_{i}$

где

${\begin{aligned}{\boldsymbol {V}}_{i}&={\begin{bmatrix}{\boldsymbol {v}}_{1}&{\boldsymbol {v}}_{2}&\cdots &{\boldsymbol {v}}_{i}\end{bmatrix}}{\text{,}}\\{\boldsymbol {\tilde {H}}}_{i}&={\begin{bmatrix}h_{11}&h_{12}&h_{13}&\cdots &h_{1,i}\\h_{21}&h_{22}&h_{23}&\cdots &h_{2,i}\\&h_{32}&h_{33}&\cdots &h_{3,i}\\&&\ddots &\ddots &\vdots \\&&&h_{i,i-1}&h_{i,i}\\&&&&h_{i+1,i}\end{bmatrix}}={\begin{bmatrix}{\boldsymbol {H}}_{i}\\h_{i+1,i}{\boldsymbol {e}}_{i}^{\mathrm {T} }\end{bmatrix}}\end{aligned}}$

с

$h_{ji}={\begin{cases}{\boldsymbol {v}}_{j}^{\mathrm {T} }{\boldsymbol {Av}}_{i}&{\text{if }}j\leq i{\text{,}}\\\lVert {\boldsymbol {w}}_{i+1}\rVert _{2}&{\text{if }}j=i+1{\text{,}}\\0&{\text{if }}j>i+1{\text{.}}\end{cases}}$

Применяя итерацию Арнольди к решению линейных систем, мы начинаем с ${\boldsymbol {r}}_{0}={\boldsymbol {b}}-{\boldsymbol {Ax}}_{0}$, остаток, соответствующий первоначальному предположению ${\boldsymbol {x}}_{0}$. После каждого шага итерации вычисляется ${\boldsymbol {y}}_{i}={\boldsymbol {H}}_{i}^{-1}(\lVert {\boldsymbol {r}}_{0}\rVert _{2}{\boldsymbol {e}}_{1})$ и новая итерация ${\boldsymbol {x}}_{i}={\boldsymbol {x}}_{0}+{\boldsymbol {V}}_{i}{\boldsymbol {y}}_{i}$.

Прямой метод Ланцоша

В ходе дальнейшего обсуждения мы предполагаем, что ${\boldsymbol {A}}$ является симметричным положительно определенным. С симметрией ${\boldsymbol {A}}$, верхняя матрица Хессенберга ${\boldsymbol {H}}_{i}={\boldsymbol {V}}_{i}^{\mathrm {T} }{\boldsymbol {AV}}_{i}$ становится симметричным и, следовательно, трехдиагональным. Тогда это можно более четко обозначить через

${\boldsymbol {H}}_{i}={\begin{bmatrix}a_{1}&b_{2}\\b_{2}&a_{2}&b_{3}\\&\ddots &\ddots &\ddots \\&&b_{i-1}&a_{i-1}&b_{i}\\&&&b_{i}&a_{i}\end{bmatrix}}{\text{.}}$

Это обеспечивает короткий трехкратный рецидив для ${\boldsymbol {v}}_{i}$ в итерации, а итерация Арнольди сводится к итерации Ланцоша.

С ${\boldsymbol {A}}$ является симметричным положительно определенным, поэтому ${\boldsymbol {H}}_{i}$. Следовательно, ${\boldsymbol {H}}_{i}$ может быть LU-факторизован без частичного поворота в

${\boldsymbol {H}}_{i}={\boldsymbol {L}}_{i}{\boldsymbol {U}}_{i}={\begin{bmatrix}1\\c_{2}&1\\&\ddots &\ddots \\&&c_{i-1}&1\\&&&c_{i}&1\end{bmatrix}}{\begin{bmatrix}d_{1}&b_{2}\\&d_{2}&b_{3}\\&&\ddots &\ddots \\&&&d_{i-1}&b_{i}\\&&&&d_{i}\end{bmatrix}}$

с удобными повторениями для $c_{i}$ и $d_{i}$:

${\begin{aligned}c_{i}&=b_{i}/d_{i-1}{\text{,}}\\d_{i}&={\begin{cases}a_{1}&{\text{if }}i=1{\text{,}}\\a_{i}-c_{i}b_{i}&{\text{if }}i>1{\text{.}}\end{cases}}\end{aligned}}$

Переписать ${\boldsymbol {x}}_{i}={\boldsymbol {x}}_{0}+{\boldsymbol {V}}_{i}{\boldsymbol {y}}_{i}$ как

${\begin{aligned}{\boldsymbol {x}}_{i}&={\boldsymbol {x}}_{0}+{\boldsymbol {V}}_{i}{\boldsymbol {H}}_{i}^{-1}(\lVert {\boldsymbol {r}}_{0}\rVert _{2}{\boldsymbol {e}}_{1})\\&={\boldsymbol {x}}_{0}+{\boldsymbol {V}}_{i}{\boldsymbol {U}}_{i}^{-1}{\boldsymbol {L}}_{i}^{-1}(\lVert {\boldsymbol {r}}_{0}\rVert _{2}{\boldsymbol {e}}_{1})\\&={\boldsymbol {x}}_{0}+{\boldsymbol {P}}_{i}{\boldsymbol {z}}_{i}\end{aligned}}$

с

${\begin{aligned}{\boldsymbol {P}}_{i}&={\boldsymbol {V}}_{i}{\boldsymbol {U}}_{i}^{-1}{\text{,}}\\{\boldsymbol {z}}_{i}&={\boldsymbol {L}}_{i}^{-1}(\lVert {\boldsymbol {r}}_{0}\rVert _{2}{\boldsymbol {e}}_{1}){\text{.}}\end{aligned}}$

Теперь важно заметить, что

${\begin{aligned}{\boldsymbol {P}}_{i}&={\begin{bmatrix}{\boldsymbol {P}}_{i-1}&{\boldsymbol {p}}_{i}\end{bmatrix}}{\text{,}}\\{\boldsymbol {z}}_{i}&={\begin{bmatrix}{\boldsymbol {z}}_{i-1}\\\zeta _{i}\end{bmatrix}}{\text{.}}\end{aligned}}$

Действительно, бывают кратковременные рецидивы. ${\boldsymbol {p}}_{i}$ и $\zeta _{i}$ также:

${\begin{aligned}{\boldsymbol {p}}_{i}&={\frac {1}{d_{i}}}({\boldsymbol {v}}_{i}-b_{i}{\boldsymbol {p}}_{i-1}){\text{,}}\\\zeta _{i}&=-c_{i}\zeta _{i-1}{\text{.}}\end{aligned}}$

Используя эту формулировку, мы приходим к простой рекуррентности для ${\boldsymbol {x}}_{i}$:

${\begin{aligned}{\boldsymbol {x}}_{i}&={\boldsymbol {x}}_{0}+{\boldsymbol {P}}_{i}{\boldsymbol {z}}_{i}\\&={\boldsymbol {x}}_{0}+{\boldsymbol {P}}_{i-1}{\boldsymbol {z}}_{i-1}+\zeta _{i}{\boldsymbol {p}}_{i}\\&={\boldsymbol {x}}_{i-1}+\zeta _{i}{\boldsymbol {p}}_{i}{\text{.}}\end{aligned}}$

Приведенные выше соотношения напрямую приводят к прямому методу Ланцоша, который оказывается несколько более сложным.

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

Если мы позволим ${\boldsymbol {p}}_{i}$ для масштабирования и компенсации масштабирования постоянного коэффициента мы потенциально можем иметь более простые повторения формы:

${\begin{aligned}{\boldsymbol {x}}_{i}&={\boldsymbol {x}}_{i-1}+\alpha _{i-1}{\boldsymbol {p}}_{i-1}{\text{,}}\\{\boldsymbol {r}}_{i}&={\boldsymbol {r}}_{i-1}-\alpha _{i-1}{\boldsymbol {Ap}}_{i-1}{\text{,}}\\{\boldsymbol {p}}_{i}&={\boldsymbol {r}}_{i}+\beta _{i-1}{\boldsymbol {p}}_{i-1}{\text{.}}\end{aligned}}$

В качестве предпосылки для упрощения мы теперь выведем ортогональность ${\boldsymbol {r}}_{i}$ и сопряженность ${\boldsymbol {p}}_{i}$, то есть для $i\neq j$,

${\begin{aligned}{\boldsymbol {r}}_{i}^{\mathrm {T} }{\boldsymbol {r}}_{j}&=0{\text{,}}\\{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ap}}_{j}&=0{\text{.}}\end{aligned}}$

Остатки взаимно ортогональны, поскольку ${\boldsymbol {r}}_{i}$ по сути является кратным ${\boldsymbol {v}}_{i+1}$ поскольку для $i=0$, ${\boldsymbol {r}}_{0}=\lVert {\boldsymbol {r}}_{0}\rVert _{2}{\boldsymbol {v}}_{1}$, для $i>0$,

${\begin{aligned}{\boldsymbol {r}}_{i}&={\boldsymbol {b}}-{\boldsymbol {Ax}}_{i}\\&={\boldsymbol {b}}-{\boldsymbol {A}}({\boldsymbol {x}}_{0}+{\boldsymbol {V}}_{i}{\boldsymbol {y}}_{i})\\&={\boldsymbol {r}}_{0}-{\boldsymbol {AV}}_{i}{\boldsymbol {y}}_{i}\\&={\boldsymbol {r}}_{0}-{\boldsymbol {V}}_{i+1}{\boldsymbol {\tilde {H}}}_{i}{\boldsymbol {y}}_{i}\\&={\boldsymbol {r}}_{0}-{\boldsymbol {V}}_{i}{\boldsymbol {H}}_{i}{\boldsymbol {y}}_{i}-h_{i+1,i}({\boldsymbol {e}}_{i}^{\mathrm {T} }{\boldsymbol {y}}_{i}){\boldsymbol {v}}_{i+1}\\&=\lVert {\boldsymbol {r}}_{0}\rVert _{2}{\boldsymbol {v}}_{1}-{\boldsymbol {V}}_{i}(\lVert {\boldsymbol {r}}_{0}\rVert _{2}{\boldsymbol {e}}_{1})-h_{i+1,i}({\boldsymbol {e}}_{i}^{\mathrm {T} }{\boldsymbol {y}}_{i}){\boldsymbol {v}}_{i+1}\\&=-h_{i+1,i}({\boldsymbol {e}}_{i}^{\mathrm {T} }{\boldsymbol {y}}_{i}){\boldsymbol {v}}_{i+1}{\text{.}}\end{aligned}}$

Чтобы увидеть сопряженность ${\boldsymbol {p}}_{i}$, достаточно показать, что ${\boldsymbol {P}}_{i}^{\mathrm {T} }{\boldsymbol {AP}}_{i}$ диагональ:

${\begin{aligned}{\boldsymbol {P}}_{i}^{\mathrm {T} }{\boldsymbol {AP}}_{i}&={\boldsymbol {U}}_{i}^{-\mathrm {T} }{\boldsymbol {V}}_{i}^{\mathrm {T} }{\boldsymbol {AV}}_{i}{\boldsymbol {U}}_{i}^{-1}\\&={\boldsymbol {U}}_{i}^{-\mathrm {T} }{\boldsymbol {H}}_{i}{\boldsymbol {U}}_{i}^{-1}\\&={\boldsymbol {U}}_{i}^{-\mathrm {T} }{\boldsymbol {L}}_{i}{\boldsymbol {U}}_{i}{\boldsymbol {U}}_{i}^{-1}\\&={\boldsymbol {U}}_{i}^{-\mathrm {T} }{\boldsymbol {L}}_{i}\end{aligned}}$

симметричен и одновременно является нижним треугольным и, следовательно, должен быть диагональным.

Теперь мы можем вывести постоянные коэффициенты $\alpha _{i}$ и $\beta _{i}$ относительно масштабированного ${\boldsymbol {p}}_{i}$ путем исключительно навязывания ортогональности ${\boldsymbol {r}}_{i}$ и сопряженность ${\boldsymbol {p}}_{i}$.

Ввиду ортогональности ${\boldsymbol {r}}_{i}$, необходимо, чтобы ${\boldsymbol {r}}_{i+1}^{\mathrm {T} }{\boldsymbol {r}}_{i}=({\boldsymbol {r}}_{i}-\alpha _{i}{\boldsymbol {Ap}}_{i})^{\mathrm {T} }{\boldsymbol {r}}_{i}=0$. Как результат,

${\begin{aligned}\alpha _{i}&={\frac {{\boldsymbol {r}}_{i}^{\mathrm {T} }{\boldsymbol {r}}_{i}}{{\boldsymbol {r}}_{i}^{\mathrm {T} }{\boldsymbol {Ap}}_{i}}}\\&={\frac {{\boldsymbol {r}}_{i}^{\mathrm {T} }{\boldsymbol {r}}_{i}}{({\boldsymbol {p}}_{i}-\beta _{i-1}{\boldsymbol {p}}_{i-1})^{\mathrm {T} }{\boldsymbol {Ap}}_{i}}}\\&={\frac {{\boldsymbol {r}}_{i}^{\mathrm {T} }{\boldsymbol {r}}_{i}}{{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ap}}_{i}}}{\text{.}}\end{aligned}}$

Аналогично, в силу сопряженности ${\boldsymbol {p}}_{i}$, необходимо, чтобы ${\boldsymbol {p}}_{i+1}^{\mathrm {T} }{\boldsymbol {Ap}}_{i}=({\boldsymbol {r}}_{i+1}+\beta _{i}{\boldsymbol {p}}_{i})^{\mathrm {T} }{\boldsymbol {Ap}}_{i}=0$. Как результат,

${\begin{aligned}\beta _{i}&=-{\frac {{\boldsymbol {r}}_{i+1}^{\mathrm {T} }{\boldsymbol {Ap}}_{i}}{{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ap}}_{i}}}\\&=-{\frac {{\boldsymbol {r}}_{i+1}^{\mathrm {T} }({\boldsymbol {r}}_{i}-{\boldsymbol {r}}_{i+1})}{\alpha _{i}{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ap}}_{i}}}\\&={\frac {{\boldsymbol {r}}_{i+1}^{\mathrm {T} }{\boldsymbol {r}}_{i+1}}{{\boldsymbol {r}}_{i}^{\mathrm {T} }{\boldsymbol {r}}_{i}}}{\text{.}}\end{aligned}}$

На этом вывод завершен.

Ссылки

1. Хестенес, MR ; Штифель, Э. (декабрь 1952 г.). «Методы сопряженных градиентов для решения линейных систем» (PDF) . Журнал исследований Национального бюро стандартов . 49 (6): 409. doi : 10.6028/jres.049.044 .
2. Шевчук, Джонатан Ричард. « Введение в метод сопряженных градиентов без мучительной боли ». (1994)
3. Саад, Ю. (2003). «Глава 6: Методы подпространств Крылова, Часть I». Итерационные методы для разреженных линейных систем (2-е изд.). СИАМ. ISBN  978-0-89871-534-7 .