Метод Гаусса: система уравнений за n³

Коротко

Система из n линейных уравнений с n неизвестными решается в два прохода. Прямой ход вычитает уравнения друг из друга так, чтобы под диагональю остались нули: с помощью первого уравнения первое неизвестное убирается из всех нижних, с помощью второго - второе, и так до треугольника. Обратный ход находит неизвестные снизу вверх. Цена - тройной вложенный цикл: n³. Замер на тренажёре: 9 920 шагов при n = 20 и 541 250 при n = 80 - вчетверо больше уравнений, в 55 раз дороже. Если на диагонали встретился ноль, деление невозможно, и строки приходится менять местами.

Если ты приведёшь систему из трёх уравнений к треугольному виду, найдёшь неизвестные обратным ходом, объяснишь, откуда в цене куб, и знаешь, что делать с нулём на диагонали, - листай до «Теперь сам».

Как думаешь, во сколько раз дороже восемьдесят уравнений

Три уравнения с тремя неизвестными:

2x + 1y + 1z =  7
4x + 3y + 3z = 19
8x + 7y + 9z = 49

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

На что это похоже

Представь, что ты распутываешь клубок из нескольких ниток, и каждая переплетена со всеми. Тянуть за все сразу бесполезно. Берёшь первую нитку и аккуратно выпутываешь её из каждой следующей. Потом вторую - из оставшихся. В конце последняя нитка лежит отдельно, предпоследняя переплетена только с ней, и клубок разбирается с конца в обратном порядке.

Метод Гаусса делает то же с неизвестными: сначала каждое по очереди убирается из всех уравнений ниже, потом уравнения решаются снизу вверх, от самого простого.

Механизм: прямой ход и обратный ход

Систему записывают таблицей коэффициентов - матрицей - и столбцом правых частей; коротко пишут Ax = b, где A - матрица, x - столбец неизвестных, b - правые части. Диагональ матрицы - клетки, где номер строки равен номеру столбца. Над уравнениями разрешены два действия: из одного уравнения вычесть другое, умноженное на число, и поменять два уравнения местами. Решения от этого не меняются.

исходная система                 убираем x из 2-й и 3-й строк          убираем y из 3-й строки
2  1  1 |  7                     2  1  1 |  7                          2  1  1 |  7
4  3  3 | 19    −2 · (1)         0  1  1 |  5                          0  1  1 |  5
8  7  9 | 49    −4 · (1)         0  3  5 | 21    −3 · (2)              0  0  2 |  6

Множитель подбирается так, чтобы нужный коэффициент обнулился: во второй строке при x стоит 4, в первой 2, значит вычесть надо первую строку, умноженную на 4 / 2 = 2. Число на диагонали, на которое делят, называют ведущим элементом: здесь это 2, потом 1, потом 2.

Получилась треугольная система. Последнее уравнение содержит одно неизвестное:

2z = 6                →  z = 3
1y + 1z = 5           →  y = 5 − 3 = 2
2x + 1y + 1z = 7      →  x = (7 − 2 − 3) / 2 = 1

Это обратный ход: каждое уравнение снизу вверх даёт одно новое неизвестное, потому что все остальные в нём уже найдены.

Откуда куб. Чтобы убрать одно неизвестное из одной строки, надо пройти по всей строке - n действий. Строк, из которых его убирают, до n. И неизвестных n. Три вложенных «для каждого» - это n³; точнее, около n³ / 3 умножений, потому что строки и столбцы с каждым шагом укорачиваются. Обратный ход - два вложенных цикла, n², и на фоне прямого его не видно.

Правило

Метод Гаусса решает систему n линейных уравнений с n неизвестными за Θ(n³) арифметических действий: прямой ход для каждого столбца i вычитает из каждой нижней строки j строку i, умноженную на a[j][i] / a[i][i]; обратный ход находит неизвестные от последнего к первому за Θ(n²). Если ведущий элемент a[i][i] равен нулю, строку i меняют местами с нижней строкой, у которой в этом столбце не ноль. Если такой строки нет, у системы нет единственного решения.

Разбор: программа и замер

def gauss(a, b):
    n = len(a)
    i = 0
    while i < n:
        j = i + 1
        while j < n:
            k = a[j][i] / a[i][i]
            c = i
            while c < n:
                a[j][c] = a[j][c] - k * a[i][c]
                c = c + 1
            b[j] = b[j] - k * b[i]
            j = j + 1
        i = i + 1
    return obratnyj_hod(a, b)

def obratnyj_hod(a, b):
    n = len(a)
    x = [0.0] * n
    i = n - 1
    while i >= 0:
        s = b[i]
        c = i + 1
        while c < n:
            s = s - a[i][c] * x[c]
            c = c + 1
        x[i] = s / a[i][i]
        i = i - 1
    return x

a - список строк, b - правые части. Внутренний цикл начинается с c = i: левее в строке i уже нули, и трогать их незачем. gauss([[2, 1, 1], [4, 3, 3], [8, 7, 9]], [7, 19, 49]) возвращает [1.0, 2.0, 3.0]. Ответ дробный, хотя числа целые: деление в Python всегда даёт дробное число.

Теперь рост. В клетке (i, j) матрицы стоит (7i + 3j) mod 10 + 1, а на диагонали к этому прибавлено 10n: с такой диагональю ведущий элемент не обнуляется. Правая часть подобрана под решение из единиц:

n        шагов         во сколько раз больше
5           275
10        1 515              5,5
20        9 920              6,5
40       71 430              7,2
80      541 250              7,6

Уравнений вдвое больше - множитель подходит к восьми: это куб. Пока n мало, заметны слагаемые с n², и множитель меньше. Ответ на вопрос из начала: от двадцати к восьмидесяти - в 55 раз. При n = 1000 это миллиард действий, при n = 100 000 - 10¹⁵.

Разбор: ноль на диагонали

0x + 1y = 1
1x + 1y = 2

Решение существует и очевидно: y = 1, x = 1. Но gauss([[0, 1], [1, 1]], [1, 2]) останавливается с ошибкой «деление на ноль»: множитель для второй строки - 1 / 0. Программа сломалась не потому, что система плохая, а потому, что уравнения записаны в неудачном порядке.

Лечение - поменять уравнения местами. Порядок уравнений на решение не влияет:

1x + 1y = 2
0x + 1y = 1           теперь ведущий элемент - 1, и вычитать уже нечего

В программе перед обработкой столбца i ищут строку от i и ниже, у которой в этом столбце не ноль, и ставят её на место i. Если в столбце одни нули - неизвестное из оставшихся уравнений исчезло, и решений либо нет (осталась строка вида 0 = 5), либо бесконечно много (строка 0 = 0). Так метод Гаусса заодно отвечает на вопрос, разрешима ли система вообще.

Наивное правило «дели на диагональ» ломается на нуле. В следующем уроке окажется, что оно ломается гораздо раньше - на числе, которое не ноль, а просто маленькое.

Что ломается без этого

Первое: вычитание из уже обработанных столбцов. Цикл с c = 0 вместо c = i даёт тот же ответ и делает в полтора раза больше работы: левая часть строк и так нулевая.

Второе: обратный ход сверху вниз. Первое уравнение треугольной системы содержит все неизвестные. Начинать надо с последнего, где неизвестное одно.

Третье: изменённые входные данные. Функция выше портит a и b: после вызова в них лежит треугольная система. Если исходная матрица нужна дальше, её копируют до вызова.

Четвёртое: сравнение ведущего элемента с нулём для дробных чисел. После нескольких вычитаний на месте настоящего нуля может стоять 1e-16. Проверка == 0 его не увидит, и программа поделит на мусор. Про это - следующий урок.

Зачем это в работе

Любая задача, где много неизвестных связаны линейно, кончается системой уравнений. Токи в электрической цепи, усилия в стержнях моста, температура в узлах сетки, подгонка прямой или кривой под точки измерений - всё это Ax = b с тысячами и миллионами неизвестных.

В библиотеках это один вызов: numpy.linalg.solve(A, b). Внутри - тот же метод Гаусса с перестановками строк, написанный так, чтобы работать блоками, которые помещаются в кеш процессора. Порядок тот же, n³, а множитель меньше в сотни раз.

Куб задаёт потолок. Плотную систему на десять тысяч неизвестных решает ноутбук; на миллион - уже нет. Большие системы решаемы потому, что в их матрицах почти все числа - нули, и методы для них устроены иначе.

И ещё одна роль: симплекс-метод из раздела про линейное программирование на каждом шаге делает ровно такое исключение - вычитает строку из остальных, чтобы убрать переменную.

Теперь сам

Приведи к треугольному виду и реши:

1x + 2y = 5
3x + 4y = 11

Система из 20 уравнений стоила 9 920 шагов. Оцени цену для 160 уравнений по закону куба.
В системе из трёх уравнений после первого шага получилось: вторая строка 0 0 2 | 4, третья 0 3 1 | 5. Что делать?

Ответ. 1. Множитель для второй строки - 3 / 1 = 3. Вычитаем: (4 − 3 · 2) y = 11 − 3 · 5, то есть −2y = −4, y = 2. Обратный ход: x = 5 − 2 · 2 = 1.

Уравнений в восемь раз больше, цена - в 8³ = 512 раз: около пяти миллионов шагов. Настоящее число будет немного меньше: при n = 20 ещё заметны слагаемые с n².
Ведущий элемент второй строки - ноль. Поменять вторую и третью строки местами: получится 0 3 1 | 5 и 0 0 2 | 4 - система уже треугольная. z = 2, y = (5 − 2) / 3 = 1.

Главное

Метод Гаусса сводит систему к треугольной, вычитая уравнения друг из друга, и решает треугольную снизу вверх. Три вложенных цикла дают n³ - надёжную, но тяжёлую цену, которая и определяет, какого размера задачи решаются напрямую. Единственное место, где алгоритм может остановиться, - деление на ведущий элемент, и перестановка строк превращает эту слабость в проверку: если переставлять нечего, у системы нет единственного решения.