Математика / Численное решение линейных систем
Метод сопряжённых градиентов для системы Ax=b
Сопряжённые градиенты решают симметричную положительно определённую систему, сохраняя полезную часть прежнего направления. На каждом шаге достаточно умножения матрицы на вектор.
Формула
Для матрицы 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. Для большой системы обычно важнее достичь требуемой точности раньше этого предела. Сильно разбросанные собственные числа затрудняют этот процесс, и тогда меняют масштаб системы с помощью предобусловливания.
Как пользоваться формулой
- Выберите начальное приближение, вычислите невязку и скопируйте её в первое направление поиска.
- Получите qₖ=Apₖ; из rₖᵀrₖ и pₖᵀqₖ вычислите длину шага, затем прибавьте αₖpₖ к xₖ.
- Вычтите αₖqₖ из невязки. Если её норма мала, проверьте остановку по заново вычисленному b−Ax.
- Если нужен следующий шаг, найдите отношение квадратов норм новой и старой невязок и обновите направление.
Историческая справка
В декабре 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
Связанные формулы
Математика
Наискорейший спуск для линейной системы: точный шаг
Наискорейший спуск для симметричной положительно определённой системы выбирает длину шага по текущей невязке. Шаг минимизирует квадратичную энергию вдоль направления спуска.
Математика
Диагональное предобусловливание сопряжённых градиентов
Диагональный предобусловливатель делит компоненты невязки на диагональ матрицы. Сопряжённые градиенты после этого строят направление в масштабе, учитывающем разные диагональные коэффициенты.
Математика
Оценка ошибки решения системы через невязку и обусловленность
Относительная ошибка приближённого решения линейной системы ограничена произведением относительной невязки на число обусловленности матрицы. Оценка объясняет, почему малая невязка сама по себе не гарантирует точного решения.