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

Алгоритм математической оптимизации

Сравнение сходимости градиентного спуска с оптимальным размером шага (зеленым) и сопряженного вектора (красным) для минимизации квадратичной функции, связанной с заданной линейной системой. Сопряженный градиент, предполагая точную арифметику, сходится максимум за n шагов, где n — размер матрицы системы (здесь n  = 2).

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

Метод сопряженных градиентов также может быть использован для решения задач неограниченной оптимизации , таких как минимизация энергии . Его обычно приписывают Магнусу Хестену и Эдуарду Штифелю , ^{[1] [2],} которые запрограммировали его на Z4 , ^[3] и тщательно исследовали его. ^{[4] [5]}

Метод бисопряженного градиента обеспечивает обобщение для несимметричных матриц. Различные нелинейные методы сопряженного градиента ищут минимумы нелинейных задач оптимизации.

Описание проблемы, решаемой с помощью сопряженных градиентов

Предположим, мы хотим решить систему линейных уравнений

${\displaystyle \mathbf {A} \mathbf {x} =\mathbf {b} }$

для вектора , где известная матрица симметрична (т.е. A T ⁼ A ) , положительно определена (т.е. x ^T Ax > 0 для всех ненулевых векторов в R ⁿ ), и вещественна , и также известна. Обозначим единственное решение этой системы через . ${\displaystyle \mathbf {x} }$ ${\displaystyle n\times n}$ ${\displaystyle \mathbf {A} }$ ${\displaystyle \mathbf {x} }$ ${\displaystyle \mathbf {b} }$ ${\displaystyle \mathbf {x} _{*}}$

Вывод как прямой метод

Метод сопряженных градиентов может быть получен из нескольких различных точек зрения, включая специализацию метода сопряженных направлений для оптимизации и вариацию итерации Арнольди / Ланцоша для задач на собственные значения . Несмотря на различия в подходах, эти выводы имеют общую тему — доказательство ортогональности остатков и сопряженности направлений поиска. Эти два свойства имеют решающее значение для разработки хорошо известной краткой формулировки метода.

Мы говорим, что два ненулевых вектора u и v сопряжены (относительно ), если ${\displaystyle \mathbf {A} }$

${\displaystyle \mathbf {u} ^{\mathsf {T}}\mathbf {A} \mathbf {v} =0.}$

Так как симметрично и положительно определено, левая часть определяет скалярное произведение ${\displaystyle \mathbf {A} }$

${\displaystyle \mathbf {u} ^{\mathsf {T}}\mathbf {A} \mathbf {v} =\langle \mathbf {u} ,\mathbf {v} \rangle _{\mathbf {A} }:=\langle \mathbf {A} \mathbf {u} ,\mathbf {v} \rangle =\langle \mathbf {u} ,\mathbf {A} ^{\mathsf {T}}\mathbf {v} \rangle =\langle \mathbf {u} ,\mathbf {A} \mathbf {v} \rangle .}$

Два вектора сопряжены тогда и только тогда, когда они ортогональны относительно этого внутреннего произведения. Сопряжение является симметричным отношением: если сопряжен с , то сопряжен с . Предположим, что ${\displaystyle \mathbf {u} }$ ${\displaystyle \mathbf {v} }$ ${\displaystyle \mathbf {v} }$ ${\displaystyle \mathbf {u} }$

${\displaystyle P=\{\mathbf {p} _{1},\dots ,\mathbf {p} _{n}\}}$

есть набор взаимно сопряженных векторов относительно , ​​т.е. для всех . Тогда образует базис для , и мы можем выразить решение в этом базисе : ${\displaystyle n}$ ${\displaystyle \mathbf {A} }$ ${\displaystyle \mathbf {p} _{i}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{j}=0}$ ${\displaystyle i\neq j}$ ${\displaystyle P}$ ${\displaystyle \mathbb {R} ^{n}}$ ${\displaystyle \mathbf {x} _{*}}$ ${\displaystyle \mathbf {Ax} =\mathbf {b} }$

${\displaystyle \mathbf {x} _{*}=\sum _{i=1}^{n}\alpha _{i}\mathbf {p} _{i}\Rightarrow \mathbf {A} \mathbf {x} _{*}=\sum _{i=1}^{n}\alpha _{i}\mathbf {A} \mathbf {p} _{i}.}$

Умножение задачи слева на вектор дает ${\displaystyle \mathbf {Ax} =\mathbf {b} }$ ${\displaystyle \mathbf {p} _{k}^{\mathsf {T}}}$

${\displaystyle \mathbf {p} _{k}^{\mathsf {T}}\mathbf {b} =\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {x} _{*}=\sum _{i=1}^{n}\alpha _{i}\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{i}=\sum _{i=1}^{n}\alpha _{i}\left\langle \mathbf {p} _{k},\mathbf {p} _{i}\right\rangle _{\mathbf {A} }=\alpha _{k}\left\langle \mathbf {p} _{k},\mathbf {p} _{k}\right\rangle _{\mathbf {A} }}$

и так

${\displaystyle \alpha _{k}={\frac {\langle \mathbf {p} _{k},\mathbf {b} \rangle }{\langle \mathbf {p} _{k},\mathbf {p} _{k}\rangle _{\mathbf {A} }}}.}$

Это дает следующий метод ^[4] решения уравнения Ax = b : найти последовательность сопряженных направлений, а затем вычислить коэффициенты . ${\displaystyle n}$ ${\displaystyle \alpha _{k}}$

Как итеративный метод

Если мы тщательно выберем сопряженные векторы , то нам не понадобятся все из них для получения хорошего приближения к решению . Поэтому мы хотим рассматривать метод сопряженных градиентов как итерационный метод. Это также позволяет нам приблизительно решать системы, где n настолько велико, что прямой метод занял бы слишком много времени. ${\displaystyle \mathbf {p} _{k}}$ ${\displaystyle \mathbf {x} _{*}}$

Обозначим начальное предположение для x _∗ через x ₀ (мы можем предположить без потери общности, что x ₀ = 0 , в противном случае рассмотрим систему Az = b − Ax ₀ вместо этого). Начиная с x ₀ мы ищем решение, и на каждой итерации нам нужна метрика, чтобы сказать нам, ближе ли мы к решению x _∗ (которое нам неизвестно). Эта метрика исходит из того факта, что решение x _∗ также является единственным минимизатором следующей квадратичной функции

${\displaystyle f(\mathbf {x} )={\tfrac {1}{2}}\mathbf {x} ^{\mathsf {T}}\mathbf {A} \mathbf {x} -\mathbf {x} ^{\mathsf {T}}\mathbf {b} ,\qquad \mathbf {x} \in \mathbf {R} ^{n}\,.}$

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

${\displaystyle \mathbf {H} (f(\mathbf {x} ))=\mathbf {A} \,,}$

и то, что минимизатор (использует D f ( x )=0) решает исходную задачу, следует из его первой производной

${\displaystyle \nabla f(\mathbf {x} )=\mathbf {A} \mathbf {x} -\mathbf {b} \,.}$

Это предполагает, что первый базисный вектор p ₀ будет отрицательным градиентом f при x = x ₀ . Градиент f равен Ax − b . Начиная с начального предположения x ₀ , это означает, что мы берем p ₀ = b − Ax ₀ . Другие векторы в базисе будут сопряжены с градиентом, отсюда и название метод сопряженных градиентов . Обратите внимание, что p ₀ также является остатком, полученным на этом начальном шаге алгоритма.

Пусть r _k будет остатком на k -м шаге:

${\displaystyle \mathbf {r} _{k}=\mathbf {b} -\mathbf {Ax} _{k}.}$

Как отмечено выше, является отрицательным градиентом при , поэтому метод градиентного спуска потребует перемещения в направлении r _k . Здесь, однако, мы настаиваем на том, что направления должны быть сопряжены друг с другом. Практический способ обеспечить это — потребовать, чтобы следующее направление поиска строилось из текущего остатка и всех предыдущих направлений поиска. Ограничение сопряжения является ограничением ортонормального типа, и, следовательно, алгоритм можно рассматривать как пример ортонормализации Грама-Шмидта . Это дает следующее выражение: ${\displaystyle \mathbf {r} _{k}}$ ${\displaystyle f}$ ${\displaystyle \mathbf {x} _{k}}$ ${\displaystyle \mathbf {p} _{k}}$

${\displaystyle \mathbf {p} _{k}=\mathbf {r} _{k}-\sum _{i<k}{\frac {\mathbf {p} _{i}^{\mathsf {T}}\mathbf {A} \mathbf {r} _{k}}{\mathbf {p} _{i}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{i}}}\mathbf {p} _{i}}$

(см. рисунок в верхней части статьи, иллюстрирующий влияние ограничения сопряженности на сходимость). Следуя этому направлению, следующее оптимальное местоположение определяется как

${\displaystyle \mathbf {x} _{k+1}=\mathbf {x} _{k}+\alpha _{k}\mathbf {p} _{k}}$

с

${\displaystyle \alpha _{k}={\frac {\mathbf {p} _{k}^{\mathsf {T}}(\mathbf {b} -\mathbf {Ax} _{k})}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}}={\frac {\mathbf {p} _{k}^{\mathsf {T}}\mathbf {r} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}},}$

где последнее равенство следует из определения . Выражение для можно вывести, если подставить выражение для x _{k +1} в f и минимизировать его по ${\displaystyle \mathbf {r} _{k}}$ ${\displaystyle \alpha _{k}}$ ${\displaystyle \alpha _{k}}$

${\displaystyle {\begin{aligned}f(\mathbf {x} _{k+1})&=f(\mathbf {x} _{k}+\alpha _{k}\mathbf {p} _{k})=:g(\alpha _{k})\\g'(\alpha _{k})&{\overset {!}{=}}0\quad \Rightarrow \quad \alpha _{k}={\frac {\mathbf {p} _{k}^{\mathsf {T}}(\mathbf {b} -\mathbf {Ax} _{k})}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}}\,.\end{aligned}}}$

Полученный алгоритм

Приведенный выше алгоритм дает наиболее простое объяснение метода сопряженных градиентов. По-видимому, алгоритм в том виде, в котором он сформулирован, требует хранения всех предыдущих направлений поиска и векторов остатков, а также множества умножений матрицы на вектор, и, таким образом, может быть вычислительно затратным. Однако более детальный анализ алгоритма показывает, что ортогонален , т. е. , для i ≠ j. И -ортогонален , т. е. , для . Это можно считать, что по мере выполнения алгоритма и охватывает то же подпространство Крылова . Где образуют ортогональный базис относительно стандартного скалярного произведения, а образуют ортогональный базис относительно скалярного произведения, индуцированного . Следовательно, можно рассматривать как проекцию на подпространство Крылова. ${\displaystyle \mathbf {r} _{i}}$ ${\displaystyle \mathbf {r} _{j}}$ ${\displaystyle \mathbf {r} _{i}^{\mathsf {T}}\mathbf {r} _{j}=0}$ ${\displaystyle \mathbf {p} _{i}}$ ${\displaystyle \mathbf {A} }$ ${\displaystyle \mathbf {p} _{j}}$ ${\displaystyle \mathbf {p} _{i}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{j}=0}$ ${\displaystyle i\neq j}$ ${\displaystyle \mathbf {p} _{i}}$ ${\displaystyle \mathbf {r} _{i}}$ ${\displaystyle \mathbf {r} _{i}}$ ${\displaystyle \mathbf {p} _{i}}$ ${\displaystyle \mathbf {A} }$ ${\displaystyle \mathbf {x} _{k}}$ ${\displaystyle \mathbf {x} }$

То есть, если метод CG начинается с , то ^[6] Ниже подробно описан алгоритм решения , где — действительная, симметричная, положительно определенная матрица. Входной вектор может быть приближенным начальным решением или 0 . Это другая формулировка точной процедуры, описанной выше. ${\displaystyle \mathbf {x} _{0}=0}$ $${\displaystyle x_{k}=\mathrm {argmin} _{y\in \mathbb {R} ^{n}}{\left\{(x-y)^{\top }A(x-y):y\in \operatorname {span} \left\{b,Ab,\ldots ,A^{k-1}b\right\}\right\}}}$$ ${\displaystyle \mathbf {A} \mathbf {x} =\mathbf {b} }$ ${\displaystyle \mathbf {A} }$ ${\displaystyle \mathbf {x} _{0}}$

${\displaystyle {\begin{aligned}&\mathbf {r} _{0}:=\mathbf {b} -\mathbf {Ax} _{0}\\&{\hbox{if }}\mathbf {r} _{0}{\text{ is sufficiently small, then return }}\mathbf {x} _{0}{\text{ as the result}}\\&\mathbf {p} _{0}:=\mathbf {r} _{0}\\&k:=0\\&{\text{repeat}}\\&\qquad \alpha _{k}:={\frac {\mathbf {r} _{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 \mathbf {r} _{k+1}:=\mathbf {r} _{k}-\alpha _{k}\mathbf {Ap} _{k}\\&\qquad {\hbox{if }}\mathbf {r} _{k+1}{\text{ is sufficiently small, then exit loop}}\\&\qquad \beta _{k}:={\frac {\mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {r} _{k+1}}{\mathbf {r} _{k}^{\mathsf {T}}\mathbf {r} _{k}}}\\&\qquad \mathbf {p} _{k+1}:=\mathbf {r} _{k+1}+\beta _{k}\mathbf {p} _{k}\\&\qquad k:=k+1\\&{\text{end repeat}}\\&{\text{return }}\mathbf {x} _{k+1}{\text{ as the result}}\end{aligned}}}$

Это наиболее часто используемый алгоритм. Та же формула для β _k используется также в методе нелинейного сопряженного градиента Флетчера–Ривза .

Перезапуски

Отметим, что вычисляется методом градиентного спуска, примененным к . Установка аналогичным образом сделала бы вычисляемым методом градиентного спуска из , т. е. может использоваться как простая реализация перезапуска итераций сопряженного градиента. ^[4] Перезапуски могут замедлить сходимость, но могут улучшить устойчивость, если метод сопряженного градиента ведет себя неправильно, например, из-за ошибки округления . ${\displaystyle \mathbf {x} _{1}}$ ${\displaystyle \mathbf {x} _{0}}$ ${\displaystyle \beta _{k}=0}$ ${\displaystyle \mathbf {x} _{k+1}}$ ${\displaystyle \mathbf {x} _{k}}$

Явный расчет остатка

Формулы и , которые обе выполняются в точной арифметике, делают формулы и математически эквивалентными. Первая используется в алгоритме, чтобы избежать дополнительного умножения на , поскольку вектор уже вычислен для оценки . Последняя может быть более точной, заменяя явное вычисление неявным рекурсией, подверженной накоплению ошибок округления , и поэтому рекомендуется для периодической оценки. ^[7] ${\displaystyle \mathbf {x} _{k+1}:=\mathbf {x} _{k}+\alpha _{k}\mathbf {p} _{k}}$ ${\displaystyle \mathbf {r} _{k}:=\mathbf {b} -\mathbf {Ax} _{k}}$ ${\displaystyle \mathbf {r} _{k+1}:=\mathbf {r} _{k}-\alpha _{k}\mathbf {Ap} _{k}}$ ${\displaystyle \mathbf {r} _{k+1}:=\mathbf {b} -\mathbf {Ax} _{k+1}}$ ${\displaystyle \mathbf {A} }$ ${\displaystyle \mathbf {Ap} _{k}}$ ${\displaystyle \alpha _{k}}$ ${\displaystyle \mathbf {r} _{k+1}:=\mathbf {b} -\mathbf {Ax} _{k+1}}$

Норма остатка обычно используется для критериев остановки. Норма явного остатка обеспечивает гарантированный уровень точности как в точной арифметике, так и при наличии ошибок округления , где сходимость естественным образом застаивается. Напротив, известно, что неявный остаток продолжает уменьшаться по амплитуде значительно ниже уровня ошибок округления и, таким образом, не может использоваться для определения застоя сходимости. ${\displaystyle \mathbf {r} _{k+1}:=\mathbf {b} -\mathbf {Ax} _{k+1}}$ ${\displaystyle \mathbf {r} _{k+1}:=\mathbf {r} _{k}-\alpha _{k}\mathbf {Ap} _{k}}$

Вычисление альфы и беты

В алгоритме α _k выбирается таким образом, чтобы было ортогонально . Знаменатель упрощается из ${\displaystyle \mathbf {r} _{k+1}}$ ${\displaystyle \mathbf {r} _{k}}$

${\displaystyle \alpha _{k}={\frac {\mathbf {r} _{k}^{\mathsf {T}}\mathbf {r} _{k}}{\mathbf {r} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}}={\frac {\mathbf {r} _{k}^{\mathsf {T}}\mathbf {r} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {Ap} _{k}}}}$

поскольку . β _k выбирается таким образом, что сопряжено с . Первоначально β _k равно ${\displaystyle \mathbf {r} _{k+1}=\mathbf {p} _{k+1}-\mathbf {\beta } _{k}\mathbf {p} _{k}}$ ${\displaystyle \mathbf {p} _{k+1}}$ ${\displaystyle \mathbf {p} _{k}}$

${\displaystyle \beta _{k}=-{\frac {\mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}}}$

с использованием

${\displaystyle \mathbf {r} _{k+1}=\mathbf {r} _{k}-\alpha _{k}\mathbf {A} \mathbf {p} _{k}}$

и эквивалентно

${\displaystyle \mathbf {A} \mathbf {p} _{k}={\frac {1}{\alpha _{k}}}(\mathbf {r} _{k}-\mathbf {r} _{k+1}),}$

числитель β _k переписывается как

${\displaystyle \mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}={\frac {1}{\alpha _{k}}}\mathbf {r} _{k+1}^{\mathsf {T}}(\mathbf {r} _{k}-\mathbf {r} _{k+1})=-{\frac {1}{\alpha _{k}}}\mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {r} _{k+1}}$

потому что и ортогональны по замыслу. Знаменатель переписывается как ${\displaystyle \mathbf {r} _{k+1}}$ ${\displaystyle \mathbf {r} _{k}}$

${\displaystyle \mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}=(\mathbf {r} _{k}+\beta _{k-1}\mathbf {p} _{k-1})^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}={\frac {1}{\alpha _{k}}}\mathbf {r} _{k}^{\mathsf {T}}(\mathbf {r} _{k}-\mathbf {r} _{k+1})={\frac {1}{\alpha _{k}}}\mathbf {r} _{k}^{\mathsf {T}}\mathbf {r} _{k}}$

используя то, что направления поиска p _k сопряжены, и снова, что остатки ортогональны. Это дает β в алгоритме после отмены α _k .

Пример кода вДжулия (язык программирования)

""" сопряженный_градиент!(A, b, x)Верните решение к `A * x = b`, используя метод сопряженных градиентов."""функция conjugate_gradient! (  A :: АбстрактнаяМатрица , b :: АбстрактныйВектор , x :: АбстрактныйВектор ; tol = eps ( eltype ( b ))   ) # Инициализация остаточного вектора остаток = b - A * x       # Инициализируем вектор направления поиска search_direction = копия ( остаточная )   # Вычислить начальную квадратичную норму остатканорма ( x ) = sqrt ( сумма ( x .^ 2 ))   old_resid_norm = норма ( остаток )   # Итерировать до сходимости в то время как old_resid_norm > tol    A_направление_поиска = A * направление_поиска     step_size = old_resid_norm ^ 2 / ( направление_поиска ' * A_ направление_поиска )       # Обновить решение @. x = x + размер_шага * направление_поиска        # Обновление остатка @. остаток = остаток - размер_шага * направление_поиска        new_resid_norm = норма ( остаток )    # Обновить вектор направления поиска @. направление_поиска = остаток +      ( новая_норма_остаточного_состояния / старая_норма_остаточного_состояния ) ^ 2 * направление_поиска     # Обновить квадрат нормы остатка для следующей итерации старая_норма_оста = новая_норма_оста   конец вернуть х конец

Числовой пример

Рассмотрим линейную систему Ax = b, заданную формулой

${\displaystyle \mathbf {A} \mathbf {x} ={\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}x_{1}\\x_{2}\end{bmatrix}}={\begin{bmatrix}1\\2\end{bmatrix}},}$

мы выполним два шага метода сопряженных градиентов, начиная с первоначального предположения

${\displaystyle \mathbf {x} _{0}={\begin{bmatrix}2\\1\end{bmatrix}}}$

для того, чтобы найти приближенное решение системы.

Решение

Для справки, точное решение:

${\displaystyle \mathbf {x} ={\begin{bmatrix}{\frac {1}{11}}\\\\{\frac {7}{11}}\end{bmatrix}}\approx {\begin{bmatrix}0.0909\\\\0.6364\end{bmatrix}}}$

Наш первый шаг — вычислить остаточный вектор r _{0 ,} связанный с x ₀ . Этот остаток вычисляется по формуле r ₀ = b - Ax ₀ , и в нашем случае равен

${\displaystyle \mathbf {r} _{0}={\begin{bmatrix}1\\2\end{bmatrix}}-{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}2\\1\end{bmatrix}}={\begin{bmatrix}-8\\-3\end{bmatrix}}=\mathbf {p} _{0}.}$

Поскольку это первая итерация, мы будем использовать остаточный вектор r ₀ в качестве начального направления поиска p ₀ ; способ выбора p _k изменится в дальнейших итерациях.

Теперь вычислим скаляр α _0, используя соотношение

${\displaystyle \alpha _{0}={\frac {\mathbf {r} _{0}^{\mathsf {T}}\mathbf {r} _{0}}{\mathbf {p} _{0}^{\mathsf {T}}\mathbf {Ap} _{0}}}={\frac {{\begin{bmatrix}-8&-3\end{bmatrix}}{\begin{bmatrix}-8\\-3\end{bmatrix}}}{{\begin{bmatrix}-8&-3\end{bmatrix}}{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}-8\\-3\end{bmatrix}}}}={\frac {73}{331}}\approx 0.2205}$

Теперь мы можем вычислить x _1, используя формулу

${\displaystyle \mathbf {x} _{1}=\mathbf {x} _{0}+\alpha _{0}\mathbf {p} _{0}={\begin{bmatrix}2\\1\end{bmatrix}}+{\frac {73}{331}}{\begin{bmatrix}-8\\-3\end{bmatrix}}\approx {\begin{bmatrix}0.2356\\0.3384\end{bmatrix}}.}$

Этот результат завершает первую итерацию, результатом является «улучшенное» приближенное решение системы, x ₁ . Теперь мы можем двигаться дальше и вычислить следующий остаточный вектор r _{1 ,} используя формулу

${\displaystyle \mathbf {r} _{1}=\mathbf {r} _{0}-\alpha _{0}\mathbf {A} \mathbf {p} _{0}={\begin{bmatrix}-8\\-3\end{bmatrix}}-{\frac {73}{331}}{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}-8\\-3\end{bmatrix}}\approx {\begin{bmatrix}-0.2810\\0.7492\end{bmatrix}}.}$

Нашим следующим шагом в этом процессе будет вычисление скаляра β ₀ , который в конечном итоге будет использоваться для определения следующего направления поиска p ₁ .

${\displaystyle \beta _{0}={\frac {\mathbf {r} _{1}^{\mathsf {T}}\mathbf {r} _{1}}{\mathbf {r} _{0}^{\mathsf {T}}\mathbf {r} _{0}}}\approx {\frac {{\begin{bmatrix}-0.2810&0.7492\end{bmatrix}}{\begin{bmatrix}-0.2810\\0.7492\end{bmatrix}}}{{\begin{bmatrix}-8&-3\end{bmatrix}}{\begin{bmatrix}-8\\-3\end{bmatrix}}}}=0.0088.}$

Теперь, используя этот скаляр β ₀ , мы можем вычислить следующее направление поиска p ₁ с помощью соотношения

${\displaystyle \mathbf {p} _{1}=\mathbf {r} _{1}+\beta _{0}\mathbf {p} _{0}\approx {\begin{bmatrix}-0.2810\\0.7492\end{bmatrix}}+0.0088{\begin{bmatrix}-8\\-3\end{bmatrix}}={\begin{bmatrix}-0.3511\\0.7229\end{bmatrix}}.}$

Теперь вычислим скаляр α _{1 ,} используя наш недавно полученный p _{1 ,} используя тот же метод, который использовался для α ₀ .

${\displaystyle \alpha _{1}={\frac {\mathbf {r} _{1}^{\mathsf {T}}\mathbf {r} _{1}}{\mathbf {p} _{1}^{\mathsf {T}}\mathbf {Ap} _{1}}}\approx {\frac {{\begin{bmatrix}-0.2810&0.7492\end{bmatrix}}{\begin{bmatrix}-0.2810\\0.7492\end{bmatrix}}}{{\begin{bmatrix}-0.3511&0.7229\end{bmatrix}}{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}-0.3511\\0.7229\end{bmatrix}}}}=0.4122.}$

Наконец, мы находим x _2, используя тот же метод, который использовался для нахождения x ₁ .

${\displaystyle \mathbf {x} _{2}=\mathbf {x} _{1}+\alpha _{1}\mathbf {p} _{1}\approx {\begin{bmatrix}0.2356\\0.3384\end{bmatrix}}+0.4122{\begin{bmatrix}-0.3511\\0.7229\end{bmatrix}}={\begin{bmatrix}0.0909\\0.6364\end{bmatrix}}.}$

Результат x ₂ является «лучшим» приближением к решению системы, чем x ₁ и x _0. Если бы в этом примере использовалась точная арифметика вместо ограниченной точности, то точное решение теоретически было бы достигнуто после n = 2 итераций ( n — порядок системы).

Свойства сходимости

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

Как итерационный метод , метод сопряженных градиентов монотонно (в энергетической норме) улучшает приближения к точному решению и может достичь требуемого допуска после относительно небольшого (по сравнению с размером задачи) числа итераций. Улучшение обычно линейно, а его скорость определяется числом обусловленности системной матрицы : чем оно больше, тем медленнее улучшение. ^[8] ${\displaystyle \mathbf {x} _{k}}$ ${\displaystyle \kappa (A)}$ ${\displaystyle A}$ ${\displaystyle \kappa (A)}$

Если велико, то обычно используется предварительная обработка для замены исходной системы на такую, которая меньше , см. ниже. ${\displaystyle \kappa (A)}$ ${\displaystyle \mathbf {Ax} -\mathbf {b} =0}$ ${\displaystyle \mathbf {M} ^{-1}(\mathbf {Ax} -\mathbf {b} )=0}$ ${\displaystyle \kappa (\mathbf {M} ^{-1}\mathbf {A} )}$ ${\displaystyle \kappa (\mathbf {A} )}$

Теорема о сходимости

Определим подмножество полиномов как

${\displaystyle \Pi _{k}^{*}:=\left\lbrace \ p\in \Pi _{k}\ :\ p(0)=1\ \right\rbrace \,,}$

где — множество многочленов максимальной степени . ${\displaystyle \Pi _{k}}$ ${\displaystyle k}$

Пусть будут итеративными приближениями точного решения , и определим ошибки как . Теперь скорость сходимости можно аппроксимировать как ^{[4] [9]} ${\displaystyle \left(\mathbf {x} _{k}\right)_{k}}$ ${\displaystyle \mathbf {x} _{*}}$ ${\displaystyle \mathbf {e} _{k}:=\mathbf {x} _{k}-\mathbf {x} _{*}}$

${\displaystyle {\begin{aligned}\left\|\mathbf {e} _{k}\right\|_{\mathbf {A} }&=\min _{p\in \Pi _{k}^{*}}\left\|p(\mathbf {A} )\mathbf {e} _{0}\right\|_{\mathbf {A} }\\&\leq \min _{p\in \Pi _{k}^{*}}\,\max _{\lambda \in \sigma (\mathbf {A} )}|p(\lambda )|\ \left\|\mathbf {e} _{0}\right\|_{\mathbf {A} }\\&\leq 2\left({\frac {{\sqrt {\kappa (\mathbf {A} )}}-1}{{\sqrt {\kappa (\mathbf {A} )}}+1}}\right)^{k}\ \left\|\mathbf {e} _{0}\right\|_{\mathbf {A} }\\&\leq 2\exp \left({\frac {-2k}{\sqrt {\kappa (\mathbf {A} )}}}\right)\ \left\|\mathbf {e} _{0}\right\|_{\mathbf {A} }\,,\end{aligned}}}$

где обозначает спектр , а обозначает число обусловленности . ${\displaystyle \sigma (\mathbf {A} )}$ ${\displaystyle \kappa (\mathbf {A} )}$

Это показывает, что итераций достаточно, чтобы уменьшить ошибку до для любого . ${\displaystyle k={\tfrac {1}{2}}{\sqrt {\kappa (\mathbf {A} )}}\log \left(\left\|\mathbf {e} _{0}\right\|_{\mathbf {A} }\varepsilon ^{-1}\right)}$ ${\displaystyle 2\varepsilon }$ ${\displaystyle \varepsilon >0}$

Обратите внимание, важный предел, когда имеет тенденцию к ${\displaystyle \kappa (\mathbf {A} )}$ ${\displaystyle \infty }$

${\displaystyle {\frac {{\sqrt {\kappa (\mathbf {A} )}}-1}{{\sqrt {\kappa (\mathbf {A} )}}+1}}\approx 1-{\frac {2}{\sqrt {\kappa (\mathbf {A} )}}}\quad {\text{for}}\quad \kappa (\mathbf {A} )\gg 1\,.}$

Этот предел показывает более высокую скорость сходимости по сравнению с итеративными методами Якоби или Гаусса–Зейделя, которые масштабируются как . ${\displaystyle \approx 1-{\frac {2}{\kappa (\mathbf {A} )}}}$

В теореме о сходимости не предполагается никакой ошибки округления , но граница сходимости обычно справедлива на практике, как теоретически объяснено ^[5] Энн Гринбаум .

Практическая конвергенция

Если инициализировано случайным образом, первый этап итераций часто является самым быстрым, поскольку ошибка устраняется в подпространстве Крылова, которое изначально отражает меньшее эффективное число условий. Второй этап сходимости обычно хорошо определяется теоретической сходимостью, связанной с , но может быть сверхлинейным, в зависимости от распределения спектра матрицы и спектрального распределения ошибки. ^[5] На последнем этапе достигается наименьшая достижимая точность, и сходимость останавливается или метод может даже начать расходиться. В типичных научных вычислительных приложениях в формате двойной точности с плавающей точкой для матриц больших размеров метод сопряженных градиентов использует критерий останова с допуском, который прекращает итерации на первом или втором этапе. ${\textstyle {\sqrt {\kappa (\mathbf {A} )}}}$ ${\displaystyle A}$

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

В большинстве случаев предобуславливание необходимо для обеспечения быстрой сходимости метода сопряженных градиентов. Если симметрично, положительно определено и имеет лучшее число обусловленности, чем , можно использовать предобуславливающий метод сопряженных градиентов. Он имеет следующий вид: ^[10] ${\displaystyle \mathbf {M} ^{-1}}$ ${\displaystyle \mathbf {M} ^{-1}\mathbf {A} }$ ${\displaystyle \mathbf {A} }$

${\displaystyle \mathbf {r} _{0}:=\mathbf {b} -\mathbf {Ax} _{0}}$

${\displaystyle {\textrm {Solve:}}\mathbf {M} \mathbf {z} _{0}:=\mathbf {r} _{0}}$

${\displaystyle \mathbf {p} _{0}:=\mathbf {z} _{0}}$

${\displaystyle k:=0\,}$

повторить

${\displaystyle \alpha _{k}:={\frac {\mathbf {r} _{k}^{\mathsf {T}}\mathbf {z} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {Ap} _{k}}}}$

${\displaystyle \mathbf {x} _{k+1}:=\mathbf {x} _{k}+\alpha _{k}\mathbf {p} _{k}}$

${\displaystyle \mathbf {r} _{k+1}:=\mathbf {r} _{k}-\alpha _{k}\mathbf {Ap} _{k}}$

если r _{k +1} достаточно мало , то выходим из цикла конец, если

${\displaystyle \mathrm {Solve} \ \mathbf {M} \mathbf {z} _{k+1}:=\mathbf {r} _{k+1}}$

${\displaystyle \beta _{k}:={\frac {\mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {z} _{k+1}}{\mathbf {r} _{k}^{\mathsf {T}}\mathbf {z} _{k}}}}$

${\displaystyle \mathbf {p} _{k+1}:=\mathbf {z} _{k+1}+\beta _{k}\mathbf {p} _{k}}$

${\displaystyle k:=k+1\,}$

конец повтора

Результат: x _{k +1}

Приведенная выше формулировка эквивалентна применению метода регулярного сопряженного градиента к предварительно обусловленной системе ^[11]

${\displaystyle \mathbf {E} ^{-1}\mathbf {A} (\mathbf {E} ^{-1})^{\mathsf {T}}\mathbf {\hat {x}} =\mathbf {E} ^{-1}\mathbf {b} }$

где

${\displaystyle \mathbf {EE} ^{\mathsf {T}}=\mathbf {M} ,\qquad \mathbf {\hat {x}} =\mathbf {E} ^{\mathsf {T}}\mathbf {x} .}$

Разложение Холецкого предобуславливателя должно использоваться для сохранения симметрии (и положительной определенности) системы. Однако это разложение не нужно вычислять, и достаточно знать . Можно показать, что имеет тот же спектр, что и . ${\displaystyle \mathbf {M} ^{-1}}$ ${\displaystyle \mathbf {E} ^{-1}\mathbf {A} (\mathbf {E} ^{-1})^{\mathsf {T}}}$ ${\displaystyle \mathbf {M} ^{-1}\mathbf {A} }$

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

Примером часто используемого предобуславливателя является неполная факторизация Холецкого . ^[12]

Использование прекондиционера на практике

Важно помнить, что мы не хотим явно инвертировать матрицу, чтобы использовать ее в процессе, поскольку инвертирование заняло бы больше времени/вычислительных ресурсов, чем решение самого алгоритма сопряженного градиента. В качестве примера предположим, что мы используем предобуславливатель, полученный из неполной факторизации Холецкого. Результирующая матрица — нижняя треугольная матрица , а матрица предобуславливателя: ${\displaystyle \mathbf {M} }$ ${\displaystyle \mathbf {M} ^{-1}}$ ${\displaystyle \mathbf {M} }$ ${\displaystyle \mathbf {L} }$

${\displaystyle \mathbf {M} =\mathbf {LL} ^{\mathsf {T}}}$

Тогда нам нужно решить:

${\displaystyle \mathbf {Mz} =\mathbf {r} }$

${\displaystyle \mathbf {z} =\mathbf {M} ^{-1}\mathbf {r} }$

Но:

${\displaystyle \mathbf {M} ^{-1}=(\mathbf {L} ^{-1})^{\mathsf {T}}\mathbf {L} ^{-1}}$

Затем:

${\displaystyle \mathbf {z} =(\mathbf {L} ^{-1})^{\mathsf {T}}\mathbf {L} ^{-1}\mathbf {r} }$

Возьмем промежуточный вектор : ${\displaystyle \mathbf {a} }$

${\displaystyle \mathbf {a} =\mathbf {L} ^{-1}\mathbf {r} }$

${\displaystyle \mathbf {r} =\mathbf {L} \mathbf {a} }$

Поскольку и и известны, и является нижним треугольным, решение для легко и вычислительно дешево с помощью прямой подстановки . Затем мы подставляем в исходное уравнение: ${\displaystyle \mathbf {r} }$ ${\displaystyle \mathbf {L} }$ ${\displaystyle \mathbf {L} }$ ${\displaystyle \mathbf {a} }$ ${\displaystyle \mathbf {a} }$

${\displaystyle \mathbf {z} =(\mathbf {L} ^{-1})^{\mathsf {T}}\mathbf {a} }$

${\displaystyle \mathbf {a} =\mathbf {L} ^{\mathsf {T}}\mathbf {z} }$

Поскольку и известны, а является верхнетреугольным, решение для легко и недорого вычисляется с помощью обратной подстановки . ${\displaystyle \mathbf {a} }$ ${\displaystyle \mathbf {L} ^{\mathsf {T}}}$ ${\displaystyle \mathbf {L} ^{\mathsf {T}}}$ ${\displaystyle \mathbf {z} }$

Используя этот метод, нет необходимости в инвертировании или явном указании чего-либо, и мы все равно получаем . ${\displaystyle \mathbf {M} }$ ${\displaystyle \mathbf {L} }$ ${\displaystyle \mathbf {z} }$

Гибкий метод предобусловленных сопряженных градиентов

В численно сложных приложениях используются сложные предобуславливатели, которые могут привести к переменной предобуславливающей, изменяющейся между итерациями. Даже если предобуславливатель симметричен и положительно определен на каждой итерации, тот факт, что он может меняться, делает приведенные выше аргументы недействительными, а в практических тестах приводит к значительному замедлению сходимости представленного выше алгоритма. Использование формулы Полака–Рибьера

${\displaystyle \beta _{k}:={\frac {\mathbf {r} _{k+1}^{\mathsf {T}}\left(\mathbf {z} _{k+1}-\mathbf {z} _{k}\right)}{\mathbf {r} _{k}^{\mathsf {T}}\mathbf {z} _{k}}}}$

вместо формулы Флетчера–Ривза

${\displaystyle \beta _{k}:={\frac {\mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {z} _{k+1}}{\mathbf {r} _{k}^{\mathsf {T}}\mathbf {z} _{k}}}}$

может значительно улучшить сходимость в этом случае. ^[13] Эту версию метода предобусловленного сопряженного градиента можно назвать ^[14] гибкой , поскольку она допускает переменную предобуславливание. Также показано ^[15], что гибкая версия является надежной, даже если предобуславливатель не является симметричным положительно определенным (SPD).

Реализация гибкой версии требует хранения дополнительного вектора. Для фиксированного предобуславливателя SPD обе формулы для β _k эквивалентны в точной арифметике, т.е. без ошибки округления . ${\displaystyle \mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {z} _{k}=0,}$

Математическое объяснение лучшей сходимости метода с формулой Полака–Рибьера заключается в том, что метод локально оптимален в этом случае, в частности, он не сходится медленнее, чем локально оптимальный метод наискорейшего спуска. ^[16]

По сравнению с локально оптимальным методом наискорейшего спуска

В обоих исходных и предобусловленных методах сопряженных градиентов нужно только задать , чтобы сделать их локально оптимальными, используя линейный поиск , методы наискорейшего спуска . При такой замене векторы p всегда совпадают с векторами z , поэтому нет необходимости хранить векторы p . Таким образом, каждая итерация этих методов наискорейшего спуска немного дешевле по сравнению с методами сопряженных градиентов. Однако последние сходятся быстрее, если только не используется (сильно) переменный и/или не-SPD предобуславливатель , см. выше. ${\displaystyle \beta _{k}:=0}$

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

Метод сопряженных градиентов также может быть получен с использованием теории оптимального управления . ^[17] При таком подходе метод сопряженных градиентов выступает в качестве оптимального регулятора обратной связи для системы с двойным интегратором . Величины и являются переменными коэффициентами усиления обратной связи. ^[17] $${\displaystyle u=k(x,v):=-\gamma _{a}\nabla f(x)-\gamma _{b}v}$$ $${\displaystyle {\dot {x}}=v,\quad {\dot {v}}=u}$$ ${\displaystyle \gamma _{a}}$ ${\displaystyle \gamma _{b}}$

Сопряженный градиент по нормальным уравнениям

Метод сопряженного градиента можно применить к произвольной матрице размером n на m , применив его к нормальным уравнениям A ^T A и вектору правой части A ^T b , поскольку A ^T A является симметричной положительно-полуопределенной матрицей для любого A . Результатом является сопряженный градиент для нормальных уравнений ( CGN или CGNR ).

А ^Т Ах = А ^Т б

Как итеративный метод, нет необходимости формировать A ^T A явно в памяти, а только выполнять умножения матрицы на вектор и транспонированной матрицы на вектор. Поэтому CGNR особенно полезен, когда A является разреженной матрицей , поскольку эти операции обычно чрезвычайно эффективны. Однако недостатком формирования нормальных уравнений является то, что число обусловленности κ( A ^T A ) равно κ ² ( A ), и поэтому скорость сходимости CGNR может быть медленной, а качество приближенного решения может быть чувствительным к ошибкам округления. Поиск хорошего предобуславливателя часто является важной частью использования метода CGNR.

Было предложено несколько алгоритмов (например, CGLS, LSQR). Алгоритм LSQR предположительно имеет лучшую численную устойчивость, когда A плохо обусловлен, т. е. A имеет большое число обусловленности .

Метод сопряженных градиентов для комплексных эрмитовых матриц

Метод сопряженных градиентов с тривиальной модификацией можно расширить до решения, при заданных комплексной матрице A и векторе b, системы линейных уравнений для комплексного вектора x, где A — эрмитова (т.е. A' = A) и положительно определенная матрица , а символ ' обозначает сопряженное транспонирование . Тривиальная модификация заключается в простой замене действительного транспонирования на сопряженное транспонирование везде. ${\displaystyle \mathbf {A} \mathbf {x} =\mathbf {b} }$

Преимущества и недостатки

Преимущества и недостатки методов сопряженных градиентов обобщены в конспектах лекций Немировски и БенТала. ^{[18] : Раздел 7.3}

Патологический пример

Этот пример взят из ^[19] Пусть , и определите Поскольку обратимо, то существует единственное решение для . Решение его методом сопряженного градиентного спуска дает нам довольно плохую сходимость: Другими словами, в процессе CG ошибка растет экспоненциально, пока внезапно не станет равной нулю, когда будет найдено единственное решение. ${\textstyle t\in (0,1)}$ $${\displaystyle W={\begin{bmatrix}t&{\sqrt {t}}&&&&\\{\sqrt {t}}&1+t&{\sqrt {t}}&&&\\&{\sqrt {t}}&1+t&{\sqrt {t}}&&\\&&{\sqrt {t}}&\ddots &\ddots &\\&&&\ddots &&\\&&&&&{\sqrt {t}}\\&&&&{\sqrt {t}}&1+t\end{bmatrix}},\quad b={\begin{bmatrix}1\\0\\\vdots \\0\end{bmatrix}}}$$ ${\displaystyle W}$ ${\textstyle Wx=b}$ $${\displaystyle \|b-Wx_{k}\|^{2}=(1/t)^{k},\quad \|b-Wx_{n}\|^{2}=0}$$

Смотрите также

• Метод бисопряженного градиента (BiCG)
• Метод сопряженных квадратов градиентов (CGS)
• Метод сопряженных остатков
• Распространение гауссовского убеждения
• Итерационный метод: линейные системы
• Подпространство Крылова
• Нелинейный метод сопряженных градиентов
• Предварительная подготовка
• Разреженное умножение матрицы на вектор

Ссылки

1. Хестенес, Магнус Р.; Штифель , Эдуард (декабрь 1952 г.). "Методы сопряженных градиентов для решения линейных систем" (PDF) . Журнал исследований Национального бюро стандартов . 49 (6): 409. doi : 10.6028/jres.049.044 .
2. Straeter, TA (1971). О расширении класса Дэвидона–Бройдена ранга один, методы минимизации квазиньютона на бесконечномерное гильбертово пространство с приложениями к задачам оптимального управления (диссертация на соискание ученой степени доктора философии). Университет штата Северная Каролина. hdl : 2060/19710026200 – через сервер технических отчетов NASA.
3. Спейзер, Амброс (2004). «Конрад Цузе и ЭРМЕТ: Ein weltweiter Architektur-Vergleich» [Конрад Цузе и ЭРМЕТ: Всемирное сравнение архитектур]. В Хеллиге, Ганс Дитер (ред.). История информатики. Visionen, Paradigmen, Leitmotive (на немецком языке). Берлин: Шпрингер. п. 185. ИСБН 3-540-00217-0.
4. Поляк, Борис (1987). Введение в оптимизацию.
5. Гринбаум, Энн (1997). Итерационные методы решения линейных систем . doi :10.1137/1.9781611970937. ISBN 978-0-89871-396-1.
6. Paquette, Elliot; Trogdon, Thomas (март 2023 г.). «Универсальность для алгоритмов сопряженного градиента и MINRES на выборочных ковариационных матрицах». Сообщения по чистой и прикладной математике . 76 (5): 1085–1136. arXiv : 2007.00640 . doi :10.1002/cpa.22081. ISSN  0010-3640.
7. Шевчук, Джонатан Р. (1994). Введение в метод сопряженных градиентов без мучительной боли (PDF) .
8. Саад, Юсеф (2003). Итерационные методы для разреженных линейных систем (2-е изд.). Филадельфия, Пенсильвания: Общество промышленной и прикладной математики. С. 195. ISBN 978-0-89871-534-7.
9. Хакбуш, В. (2016-06-21). Итеративное решение больших разреженных систем уравнений (2-е изд.). Швейцария: Springer. ISBN 978-3-319-28483-5. OCLC  952572240.
10. Барретт, Ричард; Берри, Майкл; Чан, Тони Ф.; Деммель, Джеймс; Донато, Джун; Донгарра, Джек; Эйкхаут, Виктор; Позо, Ролдан; Ромин, Чарльз; ван дер Ворст, Хенк. Шаблоны для решения линейных систем: строительные блоки для итерационных методов (PDF) (2-е изд.). Филадельфия, Пенсильвания: SIAM. стр. 13. Получено 31.03.2020 .
11. Голуб, Джин Х.; Ван Лоан, Чарльз Ф. (2013). Матричные вычисления (4-е изд.). Johns Hopkins University Press. раздел 11.5.2. ISBN 978-1-4214-0794-4.
12. Concus, P.; Golub, GH; Meurant, G. (1985). «Блочная предварительная подготовка для метода сопряженных градиентов». Журнал SIAM по научным и статистическим вычислениям . 6 (1): 220–252. doi :10.1137/0906018.
13. Голуб, Джин Х.; Йе, Цян (1999). «Неточный предобусловленный метод сопряженных градиентов с внутренней-внешней итерацией». Журнал SIAM по научным вычислениям . 21 (4): 1305. CiteSeerX 10.1.1.56.1755 . doi :10.1137/S1064827597323415. 
14. Notay, Yvan (2000). «Гибкие сопряженные градиенты». Журнал SIAM по научным вычислениям . 22 (4): 1444–1460. CiteSeerX 10.1.1.35.7473 . doi :10.1137/S1064827599362314. 
15. Bouwmeester, Henricus; Dougherty, Andrew; Knyazev, Andrew V. (2015). «Несимметричное предварительное обусловливание для методов сопряженного градиента и наискорейшего спуска 1». Procedia Computer Science . 51 : 276–285. arXiv : 1212.6680 . doi : 10.1016/j.procs.2015.05.241 . S2CID  51978658.
16. Князев, Эндрю В.; Лашук, Илья (2008). «Методы наискорейшего спуска и сопряженных градиентов с переменной предобуславливающей обработкой». Журнал SIAM по матричному анализу и приложениям . 29 (4): 1267. arXiv : math/0605767 . doi :10.1137/060675290. S2CID  17614913.
17. Росс, IM , «Оптимальная теория управления для ускоренной оптимизации», arXiv :1902.09004, 2019.
18. Немировски и Бен-Тал (2023). «Оптимизация III: Выпуклая оптимизация» (PDF) .
19. Пеннингтон, Фабиан Педрегоса, Кортни Пакетт, Том Трогдон, Джеффри. «Учебник по теории случайных матриц и машинному обучению». random-matrix-learning.github.io . Получено 05.12.2023 .{{cite web}}: CS1 maint: multiple names: authors list (link)

Дальнейшее чтение

• Аткинсон, Кенделл А. (1988). "Раздел 8.9". Введение в численный анализ (2-е изд.). John Wiley and Sons. ISBN 978-0-471-50023-0.
• Avriel, Mordecai (2003). Нелинейное программирование: анализ и методы . Dover Publishing. ISBN 978-0-486-43227-4.
• Голуб, Джин Х.; Ван Лоан, Чарльз Ф. (2013). "Глава 11". Матричные вычисления (4-е изд.). Издательство Университета Джонса Хопкинса. ISBN 978-1-4214-0794-4.
• Саад, Юсеф (2003-04-01). "Глава 6" . Итерационные методы для разреженных линейных систем (2-е изд.). SIAM. ISBN 978-0-89871-534-7.
• Жерар Мёран: «Обнаружение и исправление скрытых ошибок в алгоритме сопряженных градиентов», Numerical Algorithms, т.92 (2023), стр.869-891. url=https://doi.org/10.1007/s11075-022-01380-1
• Мёрант, Жерар; Тичи, Петр (2024). Оценка нормы ошибки в алгоритме сопряженного градиента . SIAM. ISBN 978-1-61197-785-1.

Внешние ссылки

• «Сопряженные градиенты, метод», Энциклопедия математики , EMS Press , 2001 [1994]