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

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

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

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

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

Описание проблемы, решаемой с

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

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

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

Деривация как прямой метод

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

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

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

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

С $\mathbf {A}$ симметричен и положительно определен, левая часть определяет скалярный продукт

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

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

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

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

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

Умножение проблемы слева $\mathbf {Ax} =\mathbf {b}$ с вектором $\mathbf {p} _{k}^{\mathsf {T}}$ урожайность

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

и так

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

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

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

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

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

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

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

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

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

$\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 - м шаге:

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

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

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

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

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

с

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

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

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

Получившийся алгоритм

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

То есть, если метод CG начинается с $\mathbf {x} _{0}=0$, затем

$$x_{k}=\mathrm {argmin} _{y\in \mathbb {R} ^{n}}{\left\{(x-y)^{\top }A(x-y)^{T}:y\in \operatorname {span} \left\{b,Ab,\ldots ,A^{k-1}b\right\}\right\}}$$

Алгоритм решения подробно описан ниже. $\mathbf {A} \mathbf {x} =\mathbf {b}$ где $\mathbf {A}$— действительная, симметричная, положительно определенная матрица. Входной вектор $\mathbf {x} _{0}$может быть приближенным начальным решением или 0 . Это другая формулировка точной процедуры, описанной выше.

${\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 также используется в нелинейном методе сопряженных градиентов Флетчера–Ривза .

Перезапуск

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

Явный расчет невязки

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

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

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

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

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

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

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

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

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

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

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

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

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

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

$\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 _отмены .

Пример кода на Julia (язык программирования)

""" 
 conjugate_gradient!(A, b, x) 

 Возвращает решение `A * x = b`, используя метод сопряженного градиента. 
 """ 
 функция   conjugate_gradient!   ( 
     A  ::  AbstractMatrix  ,   b  ::  AbstractVector  ,   x  ::  AbstractVector  ;   tol  =  eps  (  eltype  (  b  )) 
 ) 
     # Инициализируем вектор остатка 
     остатка   =   b   -   A   *   x 
     # Инициализируем вектор направления поиска 
     search_direction   =   остаток 
     # Вычисляем начальный квадрат остатка норма 
	 норма  (  x  )   =   sqrt  (  сумма  (  x  .^  2  )) 
     old_resid_norm   =   норма  (  остаток  ) 

     # Итерация до сходимости 
     while   old_resid_norm   >   tol 
         A_search_direction   =   A   *   search_direction 
         Step_size   =   old_resid_norm  ^  2   /   (  search_direction  '   *   A_search_direction  ) 
         # Обновить решение 
         @.    x   =   x   +   размер_шага   *   направление_поиска 
         # Обновление остатка 
         @.    Остаток   =   остаток   -   размер_шага   *   A_search_direction 
         new_resid_norm   =   норма  (  остаток  ) 
        
         # Обновить вектор направления поиска 
         @.    search_direction   =   остаток   +  
             (  new_resid_norm   /   old_resid_norm  )  ^  2   *   search_direction 
         # Обновление квадрата нормы остатка для следующей итерации 
         old_resid_norm   =   new_resid_norm 
     end 
     return   x 
 end 

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

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

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

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

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

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

Решение

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

$\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 ₀ и в нашем случае равна

$\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, используя соотношение

$\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 ₁ по формуле

$\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 ₁ по формуле

$\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 ₁ .

$\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 _1, используя соотношение

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

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

$\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 ₁ .

$\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 ₀ . Если бы в этом примере использовалась точная арифметика вместо ограниченной точности, то точное решение теоретически было бы достигнуто после n = 2 итераций ( n — порядок системы).

Свойства конвергенции

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

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

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

Теорема

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

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

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

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

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

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

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

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

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

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

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

конвергенция

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

Метод предварительно

См. также: Предварительный кондиционер

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

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

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

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

$k:=0\,$

повторить

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

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

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

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

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

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

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

$k:=k+1\,$

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

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

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

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

где

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

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

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

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

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

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

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

Тогда нам предстоит решить:

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

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

Но:

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

Затем:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

$$u=k(x,v):=-\gamma _{a}\nabla f(x)-\gamma _{b}v$$

для системы двойного интегратора ,

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

Количества $\gamma _{a}$ и $\gamma _{b}$ являются переменными коэффициентами обратной связи.

нормальных

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

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

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

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

для сложных эрмитовых

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

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

Преимущества и недостатки методов сопряженных градиентов обобщены в конспектах лекций Немировского и БенТала.

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

Этот пример из Позволять ${\textstyle t\in (0,1)}$и определить

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

С $W$ обратима, то существует единственное решение ${\textstyle Wx=b}$. Решение ее методом сопряженного градиентного спуска дает нам довольно плохую сходимость:

$$\|b-Wx_{k}\|^{2}=(1/t)^{k},\quad \|b-Wx_{n}\|^{2}=0$$

Другими словами, в процессе компьютерной графики ошибка растет экспоненциально, пока внезапно не станет нулевой, когда будет найдено уникальное решение.

См. также

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

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

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

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

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