Метод Гаусса-Зейделя СЛАУ: формулы, сходимость, сравнение с Якоби

Метод Гаусса–Зейделя - итерационный способ решения системы линейных алгебраических уравнений $Ax = b$. В отличие от прямых методов (например, исключение Гаусса с обратным ходом), он строит последовательность приближений $x^{(0)}, x^{(1)}, x^{(2)}, \ldots$, сходящуюся к точному решению, и применяется там, где матрица $A$ большая, разреженная и плохо помещается в память для прямого LU-разложения: дискретизация уравнения Лапласа, метод конечных элементов, задачи теплопроводности.

Постановка задачи

Дана система $Ax = b$ из $n$ уравнений с $n$ неизвестными, где $A \in \mathbb{R}^{n \times n}$, $b \in \mathbb{R}^n$. Считаем, что $a_{ii} \neq 0$ для всех $i = 1, \ldots, n$ - иначе соответствующее уравнение нужно переставить с другим, где диагональный элемент ненулевой. Итерационные методы стартуют с произвольного начального приближения $x^{(0)}$ (обычно нулевой вектор) и на каждом шаге уточняют все компоненты вектора.

Расщепим $A$ в сумму трёх матриц: $A = L + D + U$, где $D$ - диагональ, $L$ - строго нижнетреугольная часть, $U$ - строго верхнетреугольная. Якоби решает $Dx^{(k+1)} = b - (L + U)x^{(k)}$, метод Гаусса–Зейделя - $(D + L)x^{(k+1)} = b - Ux^{(k)}$. Разница в одном символе перевешивается тем, что Гаусс–Зейдель использует уже обновлённые компоненты текущей итерации.

Идея последовательного обновления

Покомпонентная формула метода Гаусса–Зейделя:

$$x_i^{(k+1)} = \frac{1}{a_{ii}}\left(b_i - \sum_{j<i} a_{ij} x_j^{(k+1)} - \sum_{j>i} a_{ij} x_j^{(k)}\right), \quad i = 1, 2, \ldots, n.$$

Ключевой момент - в первой сумме под $j < i$ стоят уже пересчитанные на текущем шаге компоненты $x_j^{(k+1)}$. Когда мы доходим до $x_2^{(k+1)}$, у нас уже есть $x_1^{(k+1)}$, и мы подставляем именно его, а не «старое» значение $x_1^{(k)}$. На $x_3^{(k+1)}$ доступны $x_1^{(k+1)}$ и $x_2^{(k+1)}$, и так далее. По сути, мы переписываем массив на месте - это удобно с точки зрения памяти: достаточно одного вектора $x$ длины $n$, тогда как Якоби требует двух (старого и нового).

Сравнение с методом Якоби

Метод Якоби в покомпонентной форме:

$$x_i^{(k+1)} = \frac{1}{a_{ii}}\left(b_i - \sum_{j \neq i} a_{ij} x_j^{(k)}\right).$$

Все компоненты считаются на основе одного и того же «старого» вектора $x^{(k)}$, поэтому порядок обновления неважен и итерация легко распараллеливается. Однако «свежая» информация о $x_j^{(k+1)}$, уже посчитанной на этом же шаге, игнорируется. Гаусс–Зейдель эту информацию использует и, как правило, сходится примерно вдвое быстрее по числу итераций для одного и того же класса матриц. Платой становится последовательность вычислений: при простой реализации компоненты считаются строго по порядку $i = 1, 2, \ldots, n$, и параллелить «в лоб» нельзя. На практике используют red–black упорядочение (как в шахматной доске), чтобы внутри одного «цвета» компоненты были независимы и считались параллельно.

Достаточные условия сходимости

Метод Гаусса–Зейделя сходится не для любой невырожденной матрицы. Два часто проверяемых достаточных условия:

• Диагональное преобладание по строкам: $|a_{ii}| > \sum_{j \neq i} |a_{ij}|$ для всех $i$. Если выполнено строго хотя бы для одной строки и нестрого для остальных, метод сходится при любом стартовом $x^{(0)}$.
• Симметрия и положительная определённость $A$ (то есть $A = A^\top$ и $x^\top A x > 0$ для всех ненулевых $x$). В этом случае Гаусс–Зейдель сходится монотонно по энергетической норме, чего нельзя гарантировать для Якоби.

Если ни одно из условий не выполнено, нужно смотреть спектральный радиус матрицы перехода. Записав одну итерацию как $x^{(k+1)} = T x^{(k)} + c$, для Гаусса–Зейделя получим матрицу перехода $T_{GS} = -(D + L)^{-1} U$. Метод сходится тогда и только тогда, когда $\rho(T_{GS}) < 1$, где $\rho$ - спектральный радиус (модуль наибольшего по модулю собственного значения). Чем меньше $\rho$, тем быстрее сходимость: число итераций для уменьшения ошибки в $10$ раз приблизительно равно $1 / \log_{10}(1/\rho)$.

Скорость сходимости и спектральный радиус

Для метода Якоби матрица перехода - $T_J = -D^{-1}(L + U)$, для Гаусса–Зейделя - $T_{GS} = -(D + L)^{-1} U$. Для широкого класса задач (так называемые согласованно упорядоченные матрицы - например, дискретизация Лапласа на регулярной сетке) выполняется красивое соотношение Янга: $\rho(T_{GS}) = \rho(T_J)^2$. Это значит, что одна итерация Гаусса–Зейделя даёт такое же уменьшение ошибки, как две итерации Якоби. Для матрицы $100 \times 100$, типичной для конечно-разностной дискретизации, это разница в десятки тысяч итераций.

Но даже Гаусс–Зейдель плох, когда $\rho(T_{GS})$ близко к $1$. Для дискретизации уравнения Пуассона на сетке $N \times N$ имеем $\rho \approx 1 - O(h^2)$, где $h = 1/N$ - шаг сетки. Удваивая разрешение, мы вчетверо увеличиваем число итераций для достижения той же точности - это уже неприемлемо для больших сеток, и нужны методы посерьёзнее: SOR, сопряжённые градиенты, многосеточные методы.

SOR-релаксация и параметр $\omega$

Метод последовательной верхней релаксации (Successive Over-Relaxation, SOR) - естественное обобщение Гаусса–Зейделя с одним параметром $\omega$. Формула:

$$x_i^{(k+1)} = (1 - \omega) x_i^{(k)} + \frac{\omega}{a_{ii}}\left(b_i - \sum_{j<i} a_{ij} x_j^{(k+1)} - \sum_{j>i} a_{ij} x_j^{(k)}\right).$$

При $\omega = 1$ это в точности Гаусс–Зейдель. При $\omega > 1$ (over-relaxation) мы «перешагиваем» через точку, которую дал бы Гаусс–Зейдель, в надежде быстрее догнать предел. При $\omega < 1$ (under-relaxation) - наоборот, делаем меньшие шаги; иногда это спасает от расходимости. Метод сходится при $0 < \omega < 2$ (необходимое условие - Острогорский–Райх), а для симметричных положительно определённых $A$ - это же условие и достаточное.

Оптимальный $\omega$ для модельных задач

Для согласованно упорядоченных матриц оптимальное значение даётся формулой:

$$\omega_{opt} = \frac{2}{1 + \sqrt{1 - \rho(T_J)^2}}.$$

При оптимальном $\omega_{opt}$ спектральный радиус SOR равен $\rho_{SOR} = \omega_{opt} - 1$. Для дискретизации Пуассона на сетке $N \times N$ получаем $\omega_{opt} \approx 2 - 2\pi/N$, и число итераций для нужной точности падает с $O(N^2)$ (Гаусс–Зейдель) до $O(N)$ - на сетке $1000 \times 1000$ это разница между миллионом и тысячей итераций.

На практике $\rho(T_J)$ часто неизвестен, и $\omega$ подбирают численно: пробуют несколько значений в диапазоне $1.0$–$1.95$ с шагом $0.05$ и смотрят, на каком метод сходится за меньшее число итераций.

Типовые применения

Гаусс–Зейдель и его варианты популярны в следующих задачах:

• Метод конечных элементов (МКЭ) для эллиптических уравнений (упругость, теплопроводность, диффузия). Матрица жёсткости разреженная, симметричная положительно определённая.
• Метод конечных разностей для уравнения Лапласа / Пуассона $\nabla^2 u = f$ на регулярной сетке. Каждый узел связан с четырьмя (в 2D) или шестью (в 3D) соседями.
• Сглаживатель в многосеточных методах. Здесь Гаусс–Зейдель используется не до сходимости, а 2–4 итерации на каждом уровне сетки - он быстро гасит высокочастотные ошибки.
• Релаксация в физических симуляциях: распределение электрического потенциала, гидродинамика несжимаемой жидкости (через уравнение для давления).

Типовые учебные задачи

В курсах численных методов чаще всего встречаются три задачи: вручную провести 2–3 итерации Гаусса–Зейделя на системе $3 \times 3$ или $4 \times 4$, проверить выполнение условия диагонального преобладания и сравнить число итераций до достижения точности $10^{-4}$ с методом Якоби на той же системе. Полезно также проследить, как меняется $\rho(T_{GS})$ для одной и той же системы при перестановке строк - диагональное преобладание зависит от порядка уравнений.

Частые ошибки

• Применить метод к системе с нулём на диагонали - без предварительной перестановки строк формула $1/a_{ii}$ даст деление на ноль.
• Перепутать, какие компоненты «свежие», и считать всё с предыдущего вектора - это превратит метод обратно в Якоби.
• Объявить метод расходящимся после 5–10 итераций. Для матриц со спектральным радиусом $0.95$ нужны сотни итераций даже для умеренной точности.
• Использовать SOR с $\omega \geq 2$ - метод гарантированно расходится.
• Брать критерий остановки $|x^{(k+1)} - x^{(k)}| < \varepsilon$ без нормировки. На больших $|x|$ относительная ошибка окажется заметной, нужно $\|x^{(k+1)} - x^{(k)}\| / \|x^{(k+1)}\| < \varepsilon$.

FAQ

Почему Гаусс–Зейдель быстрее Якоби, если формула почти та же? Каждая компонента $x_i^{(k+1)}$ считается уже с использованием обновлённых $x_1^{(k+1)}, \ldots, x_{i-1}^{(k+1)}$. Это «свежая» информация, и она ускоряет распространение поправок по системе. Для модельных задач спектральный радиус матрицы перехода в квадрате меньше, что эквивалентно удвоенной скорости сходимости.

Что делать, если диагональное преобладание не выполнено? Сначала попробовать переставить строки (или строки и столбцы) так, чтобы максимальные по модулю элементы оказались на диагонали - это аналог поиска ведущего элемента. Если и после этого преобладания нет, проверить, не симметрична ли матрица и не положительно ли определена - этого тоже достаточно. В крайнем случае - посчитать $\rho(T_{GS})$ численно и убедиться, что он меньше $1$.

Можно ли распараллелить Гаусс–Зейдель? Не «в лоб» - компоненты считаются последовательно. Но можно использовать red–black упорядочение: в дискретизации Лапласа разбить узлы на «красные» и «чёрные» в шахматном порядке. Внутри одного цвета узлы между собой не связаны, и пересчёт всех «красных» - параллельный, потом всех «чёрных» - тоже параллельный. На GPU это даёт значительное ускорение.

Коротко

Метод Гаусса–Зейделя - итерационный метод для $Ax = b$, использующий уже обновлённые компоненты на текущем шаге. Достаточные условия сходимости - диагональное преобладание или симметричная положительно определённая матрица. Сходится примерно вдвое быстрее Якоби по числу итераций. Обобщение с параметром $\omega$ (SOR) при оптимальном выборе радикально ускоряет метод для задач с большим спектральным радиусом, типичных для дискретизации Лапласа.