Математика / Численное решение линейных систем
GMRES(1): шаг с минимальной евклидовой невязкой
GMRES с перезапуском после одного шага выбирает поправку вдоль текущей невязки так, чтобы минимизировать норму нового остатка. Формулу можно применять и к несимметричной матрице.
Формула
У кривой для примера минимум 1/5 при s=2/5. При повороте минимум остаётся в начальной точке.
Обозначения
- $A$
- вещественная квадратная матрица, —
- $b$
- правая часть системы, —
- $xₖ$
- приближение на начале цикла, —
- $rₖ$
- невязка b−Axₖ, —
- $k$
- номер одношагового цикла, —
Условия применения
На каждом цикле допустимы точки xₖ+s rₖ с вещественным s. Невязка rₖ ненулевая и Arₖ≠0; второе автоматически выполняется для обратимой A. Если rₖ=0, решение найдено, а следующий цикл выполнять не надо.
Ограничения
Минимизация даёт невозрастание нормы, но не гарантирует строгое уменьшение: при (Arₖ)ᵀrₖ=0 шаг нулевой. В полном GMRES или при перезапуске через m>1 шагов пространство поиска больше; формулу нельзя повторять вместо построения этого пространства.
Подробное объяснение
Если к xₖ прибавить s rₖ, остаток изменится на −sArₖ. Его можно понимать как вектор rₖ, из которого пытаются убрать часть, направленную вдоль q=Arₖ. Лучшая поправка в обычной евклидовой норме — ортогональная проекция rₖ на прямую q.
Квадрат нормы rₖ−sq равен rₖᵀrₖ−2s qᵀrₖ+s² qᵀq. Коэффициент при квадрате положителен, поэтому решение уравнения производной даёт минимум. После шага новый остаток ортогонален q; ортогональность к старой невязке здесь не требуется.
В первом цикле GMRES пространство поправок натянуто на rₖ. После двух шагов без перезапуска к нему добавилось бы направление Arₖ. При перезапуске после первого шага старое пространство отбрасывается, и расчёт снова начинает строить его из очередного остатка.
Если угол между rₖ и q прямой, проекция нулевая, даже когда A обратима и решение далеко. Пример — поворот на четверть оборота в плоскости. Более широкое пространство из rₖ и Arₖ в таком случае уже охватывает плоскость и позволяет решить двумерную задачу.
Как пользоваться формулой
- Начните с прямого расчёта невязки b−Axₖ; проверьте, что она выше выбранного порога остановки.
- Умножьте невязку на A и найдите скалярное произведение полученного вектора с исходной невязкой.
- Разделите это произведение на квадрат нормы Arₖ и прибавьте к приближению полученный коэффициент, умноженный на rₖ.
- Проверьте новую невязку и её ортогональность к Arₖ; при нулевом шаге разберите возможность стагнации.
Историческая справка
Юсеф Саад и Мартин Шульц опубликовали статью о GMRES в 1986 году. Их задача состояла в распространении принципа минимального остатка на несимметричные системы. Статья формулирует поиск приближения в пространствах Крылова и обсуждает перезапуск для ограничения вычислительных затрат и памяти.
В одномерном пространстве эта постановка сводится к скалярной проекции, поэтому ортогонализацию целого базиса выполнять не приходится. Руководство Templates 1994 года отдельно обсуждает выбор длины цикла: более короткий цикл уменьшает объём хранимых данных, но отбрасывает направления, которые могли бы существенно уменьшить остаток. В предельном одношаговом случае остаётся только текущая невязка.
Историческая линия формулы
Название GMRES относится к алгоритму Саада и Шульца. Приведённый коэффициент является одномерным следствием их задачи минимизации невязки: он описывает первый шаг полного GMRES или каждый цикл версии GMRES(1). Это не формула всех последующих шагов без перезапуска.
Пример
Пусть A имеет строки (1,1) и (0,2), b=(0,1), x₀=(0,0). Найдём первый шаг минимизации невязки для этой несимметричной системы. Имеем r₀=(0,1) и q=(1,2). Скалярное произведение qᵀr₀=2, а qᵀq=5. Берём s=2/5 и получаем x₁=(0,2/5). После умножения A на x₁ остаётся r₁=(−2/5,1/5). Квадрат нормы остатка стал равен 1/5 вместо 1. Проверка оптимальности: qᵀr₁=−2/5+2/5=0. Ответ для первого шага: x₁=(0,2/5), норма невязки 1/√5. Все координаты здесь без физических единиц. Поскольку новый остаток ненулевой, систему ещё предстоит решать.
Частая ошибка
В знаменателе стоит квадрат нормы Arₖ; замена на rₖᵀArₖ меняет критерий. Не считайте отрицательный шаг ошибкой: направление на прямой допускает оба знака. Если шаг равен нулю при ненулевой невязке, проверьте угол между rₖ и Arₖ, а не повторяйте тот же расчёт без изменений.
Практика
Задачи с решением
Остаток на собственной прямой
Условие. Для A из примера, b=(2,0) и x₀=0 найдите первый шаг.
Решение. q=Ar₀=(2,0)=r₀. Числитель и знаменатель равны 4, поэтому s=1. Прямое умножение A(2,0) возвращает b.
Ответ. x₁=(2,0), остаток нулевой.
Стагнация при повороте
Условие. A поворачивает вектор (u,v) в (−v,u). При b=(1,0) и x₀=0 найдите s.
Решение. Невязка (1,0) ортогональна q=(0,1); числитель ноль, знаменатель единица. Следующий цикл стартует из той же точки.
Ответ. s=0, норма остатка остаётся равной 1.
Дополнительные источники
- https://web.stanford.edu/class/cme324/saad-schultz.pdf
- https://www.netlib.org/utk/papers/templates/report.html
Связанные формулы
Математика
Наискорейший спуск для линейной системы: точный шаг
Наискорейший спуск для симметричной положительно определённой системы выбирает длину шага по текущей невязке. Шаг минимизирует квадратичную энергию вдоль направления спуска.
Математика
Метод сопряжённых градиентов для системы Ax=b
Сопряжённые градиенты решают симметричную положительно определённую систему, сохраняя полезную часть прежнего направления. На каждом шаге достаточно умножения матрицы на вектор.
Математика
Оценка ошибки решения системы через невязку и обусловленность
Относительная ошибка приближённого решения линейной системы ограничена произведением относительной невязки на число обусловленности матрицы. Оценка объясняет, почему малая невязка сама по себе не гарантирует точного решения.