Метод сопряженных градиентов - 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} ^ {п}}$, и мы можем выразить решение Икс_∗ из ${ 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 я 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 е ( 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 набла ^ {2} е ( mathbf {x}) = mathbf {A} ,,}$

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

${ Displaystyle набла е ( 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 _ {я$

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

${ 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 {выровнен} е ( 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 {align}}}$

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

Вышеупомянутый алгоритм дает наиболее прямое объяснение метода сопряженных градиентов. По-видимому, заявленный алгоритм требует хранения всех предыдущих направлений поиска и векторов остатков, а также множества операций умножения матрицы на вектор и, следовательно, может быть дорогостоящим в вычислительном отношении. Однако более внимательный анализ алгоритма показывает, что р_я ортогонален р_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 {выровнено} & mathbf {r} _ {0}: = mathbf {b} - mathbf {Ax} _ {0} & { hbox {if}} mathbf {r } _ {0} { text {достаточно маленький, тогда верните}} mathbf {x} _ {0} { text {в качестве результата}} & 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 {достаточно мал, тогда выйдите из цикла}} & 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 {в качестве результата}} end {выравнивается}}}$

Это наиболее часто используемый алгоритм. Та же формула для β_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}} приблизительно { 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 каппа (А)}$ матрицы системы ${ displaystyle A}$: больший ${ Displaystyle каппа (А)}$ есть, тем медленнее улучшение.^[5]

Если ${ Displaystyle каппа (А)}$ большой, предварительная подготовка используется для замены оригинальной системы ${ displaystyle mathbf {Ax} - mathbf {b} = 0}$ с ${ Displaystyle mathbf {M} ^ {- 1} ( mathbf {Ax} - mathbf {b}) = 0}$ такой, что ${ Displaystyle каппа ( mathbf {M} ^ {- 1} mathbf {A})}$ меньше чем ${ Displaystyle каппа ( 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 {align} 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 { каппа ( mathbf {A})}} + 1}} right) ^ {k} left | mathbf {e} _ {0} right | _ { mathbf {A}} ,, конец {выровнено}}}$

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

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

${ displaystyle { frac {{ sqrt { kappa ( mathbf {A})}} - 1} {{ sqrt { kappa ( mathbf {A})}} + 1}} приблизительно 1- { frac {2} { sqrt { kappa ( mathbf {A})}}} quad { text {for}} quad kappa ( mathbf {A}) gg 1 ,.}$

Этот предел показывает более высокую скорость сходимости по сравнению с итерационными методами Якоби или же Гаусс-Зейдель которые масштабируются как ${ Displaystyle приблизительно 1 - { гидроразрыва {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 к: = 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 к: = к + 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 {Икс} .}$

Матрица предобуславливателя 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 и знак равно К (х, v): = - гамма _ {а} набла ф (х) - гамма _ {b} v}$$

$${ displaystyle { dot {x}} = v, quad { dot {v}} = u}$$

${ displaystyle gamma _ {a}}$

${ displaystyle gamma _ {b}}$

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

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

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

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

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

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

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

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

• Аткинсон, Кенделл А. (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]