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

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

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

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

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

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

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

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

для вектора Икс, где известные п × п матрица А является симметричный (т.е. А^Т = А), положительно определенный (т.е. Икс^ТТопор > 0 для всех ненулевых векторов Икс в р^п), и настоящий, и б также известен. Обозначим единственное решение этой системы через ${\displaystyle \mathbf {x} _{*}}$.

Как прямой метод

Мы говорим, что два ненулевых вектора ты и v находятся сопрягать (относительно А) если

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

С А симметрична и положительно определена, левая часть определяет внутренний продукт

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

Два вектора сопряжены тогда и только тогда, когда они ортогональны относительно этого внутреннего продукта. Сопряжение является симметричным отношением: если ты сопряжен с v, тогда v сопряжен с ты. Предположим, что

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

это набор п взаимно сопряженные векторы (относительно А). потом п образует основа за ${\displaystyle \mathbb {R} ^{n}}$, и мы можем выразить решение Икс_∗ из ${\displaystyle \mathbf {Ax} =\mathbf {b} }$ в этой основе:

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

На основе этого расширения мы вычисляем:

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

Умножение слева на ${\displaystyle \mathbf {p} _{k}^{\mathsf {T}}}$:

${\displaystyle \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},}$

замена ${\displaystyle \mathbf {Ax_{*}} =\mathbf {b} }$ и ${\displaystyle \mathbf {u} ^{\mathsf {T}}\mathbf {A} \mathbf {v} =\langle \mathbf {u} ,\mathbf {v} \rangle _{\mathbf {A} }}$:

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

тогда ${\displaystyle \mathbf {u} ^{\mathsf {T}}\mathbf {v} =\langle \mathbf {u} ,\mathbf {v} \rangle }$ и используя ${\displaystyle \forall i\neq k:\langle \mathbf {p} _{k},\mathbf {p} _{i}\rangle _{\mathbf {A} }=0}$ дает

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

что подразумевает

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

Это дает следующий метод решения уравнения Топор = б: найти последовательность п сопряженных направлений, а затем вычислить коэффициенты α_k.

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

Если выбрать сопряженные векторы п_k внимательно, тогда нам могут не понадобиться все они для получения хорошего приближения к решению Икс_∗. Итак, мы хотим рассматривать метод сопряженных градиентов как итерационный метод. Это также позволяет нам приближенно решать системы, в которых п настолько велик, что прямой метод занял бы слишком много времени.

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

${\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 \nabla ^{2}f(\mathbf {x} )=\mathbf {A} \,,}$

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

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

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

Позволять р_k быть остаточный на kшаг:

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

Как отмечалось выше, р_k отрицательный градиент ж в Икс = Икс_k, Итак градиентный спуск метод потребует движения в направлении р_k. Однако здесь мы настаиваем на том, чтобы направления п_k быть сопряженными друг другу. Практический способ обеспечить это - требовать, чтобы следующее направление поиска строилось из текущего остатка и всех предыдущих направлений поиска.^[4] Это дает следующее выражение:

${\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}}},}$

где последнее равенство следует из определения р_k .Выражение для ${\displaystyle \alpha _{k}}$ можно получить, если подставить выражение для Икс_k+1 в ж и минимизировать его по сравнению с ${\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}}}$

Результирующий алгоритм

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

Алгоритм подробно описан ниже для решения Топор = б куда А - вещественная симметричная положительно определенная матрица. Входной вектор Икс₀ может быть приближенным начальным решением или 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 также используется в методике Флетчера – Ривза нелинейный метод сопряженных градиентов.

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

В алгоритме α_k выбирается так, что ${\displaystyle \mathbf {r} _{k+1}}$ ортогонален р_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}}}}$

поскольку ${\displaystyle \mathbf {r} _{k+1}=\mathbf {p} _{k+1}-\mathbf {\beta } _{k}\mathbf {p} _{k}}$. В β_k выбирается так, что ${\displaystyle \mathbf {p} _{k+1}}$ сопряжен с п_k. Первоначально, β_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}}$ и р_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}}$

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

Пример кода в MATLAB / GNU Octave

функцияИкс =конград(А, б, х)р = б - А * Икс;    п = р;    продано = р' * р;    за i = 1: длина (b)        Ap = А * п;        альфа = продано / (п' * Ap);        Икс = Икс + альфа * п;        р = р - альфа * Ap;        rsnew = р' * р;        если sqrt (rsnew) <1e-10              переменаконец        п = р + (rsnew / продано) * п;        продано = rsnew;    конецконец

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

Рассмотрим линейную систему Топор = б данный

${\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}}}$

Наш первый шаг - вычислить остаточный вектор р₀ связана с Икс₀. Этот остаток вычисляется по формуле р₀ = б - Топор₀, а в нашем случае равно

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

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

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

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

Теперь мы можем вычислить Икс₁ используя формулу

${\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}}={\begin{bmatrix}0.2356\\0.3384\end{bmatrix}}.}$

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

${\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}}={\begin{bmatrix}-0.2810\\0.7492\end{bmatrix}}.}$

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

${\displaystyle \beta _{0}={\frac {\mathbf {r} _{1}^{\mathsf {T}}\mathbf {r} _{1}}{\mathbf {r} _{0}^{\mathsf {T}}\mathbf {r} _{0}}}={\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.}$

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

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

Теперь вычислим скаляр α₁ используя наши недавно приобретенные п₁ используя тот же метод, что и для α₀.

${\displaystyle \alpha _{1}={\frac {\mathbf {r} _{1}^{\mathsf {T}}\mathbf {r} _{1}}{\mathbf {p} _{1}^{\mathsf {T}}\mathbf {Ap} _{1}}}={\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.}$

Наконец, мы находим Икс₂ используя тот же метод, что и для поиска Икс₁.

${\displaystyle \mathbf {x} _{2}=\mathbf {x} _{1}+\alpha _{1}\mathbf {p} _{1}={\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}}.}$

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

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

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

Если ${\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}$.

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

${\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} }\,,\end{aligned}}}$

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

Обратите внимание на важный предел, когда ${\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} )}}}$.

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

Смотрите также: Прекондиционер

В большинстве случаев, предварительная подготовка необходимо для обеспечения быстрой сходимости метода сопряженных градиентов. Предварительно обусловленный метод сопряженных градиентов имеет следующий вид:^[7]

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

${\displaystyle \mathbf {z} _{0}:=\mathbf {M} ^{-1}\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}}$

если р_k+1 достаточно мал тогда цикл выхода конец, если

${\displaystyle \mathbf {z} _{k+1}:=\mathbf {M} ^{-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\,}$

конец повторения

Результат Икс_k+1

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

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

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

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

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

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

${\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}}}}$

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

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

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

Пример кода в MATLAB / GNU Octave

функция[x, k] =cgp(x0, A, C, b, мит, stol, bbA, bbC)% Сводка:% x0: начальная точка% A: Матрица A системы Ax = b% C: матрица предварительной подготовки может быть левой или правой% mit: максимальное количество итераций% stol: остаток нормы толерантности% bbA: Черный ящик, который вычисляет произведение матрица-вектор для A * u% bbC: Черный ящик, который вычисляет:% для левого предобуславливателя: ha = C  ra% для правого прекондиционера: ha = C * ra% x: расчетная точка решения% k: количество выполненных итераций %% Пример:% tic; [x, t] = cgp (x0, S, speye (1), b, 3000, 10 ^ -8, @ (Z, o) Z * o, @ (Z, o) o); toc% Истекшее время составляет 0,550190 секунд.%% Ссылка:% Métodos iterativos tipo Krylov para sistema lineales% Б. Молина и М. Райдан - {{ISBN | 908-261-078-X}}        если nargin <8, error ("Недостаточно входных аргументов. Попробуйте помочь."); конец;        если isempty (A), error ('Входная матрица A не должна быть пустой.'); конец;        если isempty (C), error ('Входная матрица предобуславливателя C не должна быть пустой.'); конец;        Икс = x0;        ха = 0;        л.с. = 0;        hpp = 0;        ра = 0;        rp = 0;        rpp = 0;        ты = 0;        k = 0;        ра = б - BBA(А, x0); % <--- ra = b - A * x0;        пока norm (ra, inf)> stol                ха = BBC(C, ра); % <--- ha = C  ra;                k = k + 1;                если (k == мит), предупреждение("GCP: MAXIT", «Достигнуто, конверсии нет».); возвращаться; конец;                hpp = л.с.;                rpp = rp;                л.с. = ха;                rp = ра;                т = rp' * л.с.;                если к == 1                        ты = л.с.;                ещеu = hp + (t / (rpp '* hpp)) * u;                конец;                Au = bbA (A, u); % <--- Au = A * u;                а = t / (u '* Au);                Икс = Икс + а * ты;                ра = rp - а * Au;        конец;

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

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

Вывод метода

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

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

Метод сопряженного градиента также может быть получен с использованием теория оптимального управления.^[13] В этом подходе метод сопряженных градиентов выпадает как контроллер оптимальной обратной связи,

$${\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}}$ - переменные коэффициенты усиления обратной связи.^[13]

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

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

А^ТТопор = А^Тб

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

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

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

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

Рекомендации

1. Hestenes, Магнус Р.; Штифель, Эдуард (Декабрь 1952 г.). «Методы сопряженных градиентов для решения линейных систем». Журнал исследований Национального бюро стандартов. 49 (6): 409. Дои:10.6028 / jres.049.044.
2. Straeter, T.A. (1971). «О расширении класса Дэвидона – Бройдена ранга один, методов квазиньютоновской минимизации на бесконечномерное гильбертово пространство с приложениями к задачам оптимального управления». Сервер технических отчетов НАСА. НАСА. HDL:2060/19710026200.
3. Шпайзер, Амброс (2004). "Конрад Цузе унд ди ЭРМЕТ: Ein weltweiter Architektur-Vergleich" [Конрад Цузе и ERMETH: всемирное сравнение архитектур]. В Hellige, Ганс Дитер (ред.). Geschichten der Informatik. Visionen, Paradigmen, Leitmotive (на немецком). Берлин: Springer. п. 185. ISBN  3-540-00217-0.
4. Ограничение сопряжения является ограничением ортонормированного типа, и, следовательно, алгоритм имеет сходство с Ортонормализация Грама-Шмидта.
5. Саад, Юсеф (2003). Итерационные методы для разреженных линейных систем (2-е изд.). Филадельфия, Пенсильвания: Общество промышленной и прикладной математики. стр.195. ISBN  978-0-89871-534-7.
6. Hackbusch, W. (21.06.2016). Итерационное решение больших разреженных систем уравнений (2-е изд.). Швейцария: Шпрингер. ISBN  9783319284835. OCLC  952572240.
7. Барретт, Ричард; Берри, Майкл; Чан, Тони Ф .; Деммель, Джеймс; Донато, июнь; Донгарра, Джек; Эйджхут, Виктор; Посо, Ролдан; Ромайн, Чарльз; ван дер Ворст, Хенк. Шаблоны для решения линейных систем: строительные блоки для итерационных методов (PDF) (2-е изд.). Филадельфия, Пенсильвания: SIAM. п. 13. Получено 2020-03-31.
8. Concus, P .; Голуб, Г. Х .; Меурант, Г. (1985). «Блок предварительной подготовки для метода сопряженных градиентов». Журнал SIAM по научным и статистическим вычислениям. 6 (1): 220–252. Дои:10.1137/0906018.
9. Golub, Gene H .; Е, Цян (1999). «Метод неточного предварительно обусловленного сопряженного градиента с внутренней и внешней итерацией». Журнал SIAM по научным вычислениям. 21 (4): 1305. CiteSeerX  10.1.1.56.1755. Дои:10.1137 / S1064827597323415.
10. Notay, Иван (2000). «Гибкие сопряженные градиенты». Журнал SIAM по научным вычислениям. 22 (4): 1444–1460. CiteSeerX  10.1.1.35.7473. Дои:10.1137 / S1064827599362314.
11. Хенрикус Боумистер, Эндрю Догерти, Эндрю В. Князев. Несимметричное предварительное кондиционирование для методов сопряженного градиента и наискорейшего спуска. Процедуры информатики, том 51, страницы 276-285, Elsevier, 2015. https://doi.org/10.1016/j.procs.2015.05.241
12. Князев, Андрей В .; Лашук, Илья (2008). «Наилучший спуск и методы сопряженного градиента с переменной предварительной подготовкой». Журнал SIAM по матричному анализу и приложениям. 29 (4): 1267. arXiv:математика / 0605767. Дои:10.1137/060675290. S2CID  17614913.
13. Росс, И.М., "Теория оптимального управления для ускоренной оптимизации", arXiv:1902.09004, 2019.

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

• Аткинсон, Кенделл А. (1988). «Раздел 8.9». Введение в численный анализ (2-е изд.). Джон Уайли и сыновья. ISBN  978-0-471-50023-0.
• Авриэль, Мардохей (2003). Нелинейное программирование: анализ и методы.. Dover Publishing. ISBN  978-0-486-43227-4.
• Golub, Gene H .; Ван Лоан, Чарльз Ф. (1996-10-15). «Глава 10». Матричные вычисления (3-е изд.). Издательство Университета Джона Хопкинса. ISBN  978-0-8018-5414-9.
• Саад, Юсеф (1 апреля 2003 г.). "Глава 6". Итерационные методы для разреженных линейных систем (2-е изд.). СИАМ. ISBN  978-0-89871-534-7.

внешняя ссылка

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