Математика / Численное решение линейных систем

Метод сопряжённых градиентов для системы Ax=b

Сопряжённые градиенты решают симметричную положительно определённую систему, сохраняя полезную часть прежнего направления. На каждом шаге достаточно умножения матрицы на вектор.

Опубликовано: Обновлено:

Формула

$$x_{k+1}=x_k+\alpha_kp_k,\quad \alpha_k=\frac{r_k^\top r_k}{p_k^\top A p_k},\quad p_{k+1}=r_{k+1}+\frac{r_{k+1}^\top r_{k+1}}{r_k^\top r_k}p_k$$
Схема Невязка на двух шагах
01200.511.522.5Норма невязкиШагСопр. градиентыНаискор. спуск

Для матрицы diag(1,4) второй сопряжённый шаг даёт точное решение; у спуска вдоль невязки остаётся ошибка.

Обозначения

$A$
симметричная положительно определённая матрица, —
$b$
правая часть, —
$xₖ$
приближение, —
$rₖ$
невязка b−Axₖ, —
$pₖ$
направление поиска, —
$αₖ$
длина шага вдоль pₖ, —
$k$
номер итерации, —

Условия применения

Начальная невязка r₀=b−Ax₀ задаёт p₀=r₀. Матрица вещественная симметричная и положительно определённая. При ненулевом направлении pₖ знаменатель pₖᵀApₖ положителен; нулевая невязка служит признаком остановки.

Ограничения

Для несимметричных и неопределённых матриц эти рекурсии не обеспечивают свойств метода. Завершение не более чем за n шагов относится к точной арифметике. При округлении рекурсивная невязка может отойти от b−Ax, поэтому конечный результат проверяют прямо.

Подробное объяснение

Пусть F(x)=xᵀAx/2−bᵀx. Обычный спуск минимизирует F по текущей невязке, но на следующем шаге может частично уничтожить достигнутое вдоль прежней прямой. Здесь подбирают направления, для которых pᵢᵀApⱼ=0 при i≠j. Для квадратичной формы это означает отсутствие смешанного члена.

В первом шаге p₀=r₀. После движения обновляется rₖ₊₁=rₖ−αₖApₖ. В точной арифметике невязка ортогональна уже пройденным направлениям. Поэтому в производной F вдоль pₖ числитель rₖᵀpₖ можно заменить на rₖᵀrₖ.

Для следующего направления к новой невязке прибавляют долю старого направления. Отношение квадратов норм выбрано так, чтобы pₖ₊₁ᵀApₖ обращалось в ноль. В результате второй минимум на плоскости направлений сохраняет уже выполненное условие первого минимума.

В n измерениях число ненулевых взаимно A-сопряжённых направлений не превосходит n. Для большой системы обычно важнее достичь требуемой точности раньше этого предела. Сильно разбросанные собственные числа затрудняют этот процесс, и тогда меняют масштаб системы с помощью предобусловливания.

Как пользоваться формулой

  1. Выберите начальное приближение, вычислите невязку и скопируйте её в первое направление поиска.
  2. Получите qₖ=Apₖ; из rₖᵀrₖ и pₖᵀqₖ вычислите длину шага, затем прибавьте αₖpₖ к xₖ.
  3. Вычтите αₖqₖ из невязки. Если её норма мала, проверьте остановку по заново вычисленному b−Ax.
  4. Если нужен следующий шаг, найдите отношение квадратов норм новой и старой невязок и обновите направление.

Историческая справка

В декабре 1952 года Магнус Хестенес и Эдуард Штифель опубликовали Methods of Conjugate Gradients for Solving Linear Systems в журнале Национального бюро стандартов США. В статье описан алгоритм для симметричных положительно определённых систем и его связь с другими способами построения сопряжённых направлений.

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

Историческая линия формулы

Классическую рекурсию связывают с публикацией Хестенеса и Штифеля 1952 года. Здесь дан частный вариант без предобусловливателя, с евклидовыми скалярными произведениями и начальным направлением, равным невязке. Для нелинейной оптимизации используются другие обоснования и варианты коэффициентов.

Пример

Возьмём A=diag(1,4), b=(1,2), x₀=(0,0). Система имеет две координаты, поэтому в точной арифметике достаточно не более двух шагов. В первом шаге r₀=p₀=(1,2), α₀=5/17 и x₁=(5,10)/17. Остаток r₁=(12,−6)/17. Отношение квадратов норм составляет (180/289)/5=36/289. Поэтому p₁=(240,−30)/289. Второй шаг: p₁ᵀAp₁=61200/289², откуда α₁=(180/289)/(61200/289²)=17/20. Прибавление α₁p₁ даёт x₂=(1,1/2). Проверка по исходным уравнениям: 1·1=1 и 4·(1/2)=2. Векторы заданы в абстрактных числовых координатах; итоговые координаты не имеют физических единиц.

Частая ошибка

Если в каждом шаге заменять pₖ на rₖ, получится наискорейший спуск. В знаменателе берите произведение с Apₖ, а не с Arₖ. При нулевой невязке нельзя обновлять коэффициенты через деление. Не проверяйте сходимость одним абсолютным числом без учёта масштаба правой части и точности задачи.

Практика

Задачи с решением

Проверка сопряжённости

Условие. Для p₀=(1,2) и p₁=(240,−30)/289 из примера проверьте A-сопряжённость.

Решение. У A=diag(1,4) образ Ap₁=(240,−120)/289. Умножаем его скалярно на p₀: (240−240)/289.

Ответ. p₀ᵀAp₁=0.

Матрица с одним собственным числом

Условие. Для A=3I, b=(3,−6), x₀=0 найдите α₀ и x₁.

Решение. Ap₀=3p₀, поэтому знаменатель втрое больше квадрата нормы r₀=p₀. Прибавляем треть (3,−6) к нулю.

Ответ. α₀=1/3, x₁=(1,−2).

Дополнительные источники

  • https://nvlpubs.nist.gov/nistpubs/jres/049/jresv49n6p409_A1b.pdf
  • https://www.netlib.org/utk/papers/templates/report.html
  • https://ics.uci.edu/~smyth/courses/cs274/readings/optimization/shewchuk.pdf

Связанные формулы

Математика

Наискорейший спуск для линейной системы: точный шаг

$x_{k+1}=x_k+\frac{r_k^\top r_k}{r_k^\top A r_k}\,r_k,\qquad r_k=b-Ax_k$

Наискорейший спуск для симметричной положительно определённой системы выбирает длину шага по текущей невязке. Шаг минимизирует квадратичную энергию вдоль направления спуска.

Математика

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

$Dz_k=r_k,\quad \alpha_k=\frac{r_k^\top z_k}{p_k^\top A p_k},\quad p_{k+1}=z_{k+1}+\frac{r_{k+1}^\top z_{k+1}}{r_k^\top z_k}p_k,\quad D=\operatorname{diag}(A)$

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

Математика

Оценка ошибки решения системы через невязку и обусловленность

$\frac{\|\widetilde{x}-x\|}{\|x\|}\le\kappa(A)\frac{\|r\|}{\|b\|}$

Относительная ошибка приближённого решения линейной системы ограничена произведением относительной невязки на число обусловленности матрицы. Оценка объясняет, почему малая невязка сама по себе не гарантирует точного решения.