Все формулы РФ

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

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

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

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

Ручной расчетСтраницы с задачами и решениямиЕсть историческая справкаДиагностика условий

Формула

$$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)$$
Схема Спектр после масштабирования
0245.500.511.52Собственные числа BcБольший кореньМеньший корень

Масштабирование даёт диагональ из единиц, но собственные числа расходятся при усилении связи c.

Обозначения

$D$
диагональная часть A, —
$A$
симметричная положительно определённая матрица, —
$rₖ$
текущая невязка исходной системы, —
$zₖ$
решение Dzₖ=rₖ, —
$pₖ$
направление поиска; p₀=z₀, —
$αₖ$
длина шага, —

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

Матрица A вещественная, симметричная и положительно определённая; следовательно, все элементы её диагонали положительны. D фиксируется на весь расчёт. Первое направление равно z₀, а x и r обновляются как x+αp и r−αAp.

Ограничения

Единичная диагональ не гарантирует хорошей обусловленности: сильные связи вне диагонали могут оставить очень малое собственное число. Не следует подавать D⁻¹A как обычную симметричную матрицу в реализацию CG; в исходных координатах нужна рекурсия с rᵀz.

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

Если у двух уравнений разные диагональные коэффициенты, одинаковые остатки могут требовать совсем разных поправок к неизвестным. Деление rᵢ на aᵢᵢ даёт поправку, которая занулила бы компоненту остатка при выключенных внедиагональных связях.

Математически применяется обычный CG к симметричной матрице B=D⁻¹ᐟ²AD⁻¹ᐟ² с преобразованной правой частью D⁻¹ᐟ²b. Подстановка y=D¹ᐟ²x возвращает вычисления в исходные координаты. Квадрат нормы преобразованной невязки равен rᵀD⁻¹r=rᵀz.

Поэтому в обеих формулах для коэффициентов встречается rᵀz, а умножение на матрицу по-прежнему выполняют с A. Явные матрицы D⁻¹ и D⁻¹ᐟ² не нужны: достаточно сохранить диагональ и поделить компоненты вектора.

При диагональной A симметрично масштабированная матрица равна I и решение получается за один шаг. При A с диагональю 4 и 9 и внедиагональным коэффициентом c матрица B имеет собственные числа 1±c/6. При приближении |c| к 6 меньший корень стремится к нулю, хотя диагональ B остаётся единичной.

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

  1. Сохраните положительные диагональные элементы A и вычислите r₀ из исходной, ещё не масштабированной системы.
  2. Поделите каждую компоненту r₀ на соответствующий диагональный элемент и присвойте полученный z₀ направлению p₀.
  3. Вычислите Apₖ и коэффициент шага через rₖᵀzₖ; обновите xₖ и исходную невязку rₖ.
  4. Проверьте критерий остановки по b−Ax. Для следующего шага повторите диагональное деление и пересчитайте направление.

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

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

В том же 1994 году Джонатан Шевчук изложил предобусловленные сопряжённые градиенты через преобразование исходной задачи. Такой вывод поясняет, почему матрица в преобразованных координатах симметрична, хотя произведение D⁻¹A в исходных координатах обычно несимметрично. В вычислительной записи остаётся решение с D, а квадратные корни из её элементов не требуются.

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

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

Пример

Пусть A=diag(2,8), b=(2,8) и x₀=(0,0). Вычислим первый шаг с диагональным предобусловливателем и проверим его прямой подстановкой. Здесь D=A. Невязка r₀=(2,8), следовательно z₀=p₀=(1,1). Числитель r₀ᵀz₀=2+8=10, знаменатель p₀ᵀAp₀=2+8=10. Значит α₀=1, и x₁=(1,1). Матрица A переводит найденный вектор точно в (2,8), поэтому r₁=0 и новый коэффициент направления не нужен. Абстрактная система не содержит размерных величин. Обычный первый шаг вдоль (2,8) не может дать (1,1) одним скалярным множителем.

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

Берите диагональ исходной A, а не диагональ обратной матрицы. Не заменяйте rᵀz на zᵀz: квадрат нормы возникает в координатах с половинным масштабированием. Для остановки проверяйте исходную невязку, так как малый z на строках с большой диагональю может скрыть существенный остаток.

Практика

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

Один шаг при ненулевой связи

Условие. A имеет диагональ 4 и 9, оба внедиагональных элемента равны 2. b=(4,9), x₀=0. Найдите x₁.

Решение. D⁻¹b=(1,1), A(1,1)=(6,11). Получаем α=(4+9)/(6+11)=13/17 и умножаем на p₀.

Ответ. x₁=(13/17,13/17).

Спектр после масштабирования

Условие. Для той же A найдите спектр B=D⁻¹ᐟ²AD⁻¹ᐟ² и отношение крайних чисел.

Решение. Внедиагональный элемент B равен 2/(2·3)=1/3. На направлениях (1,1) и (1,−1) получаем 1+1/3 и 1−1/3.

Ответ. 4/3 и 2/3; отношение равно 2.

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

  • https://www.netlib.org/utk/papers/templates/report.html
  • https://ics.uci.edu/~smyth/courses/cs274/readings/optimization/shewchuk.pdf

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

Математика

Метод сопряжённых градиентов для системы 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$

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

Математика

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

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

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

Математика

Метод Якоби для системы линейных уравнений

$x_i^{(k+1)}=\frac{b_i-\sum_{j\ne i}a_{ij}x_j^{(k)}}{a_{ii}},\quad i=1,\ldots,n$

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