NUM Числа и алгебра Алгоритм
Метод Гаусса
Также: Gaussian elimination · LU decomposition · LU-разложение · метод исключения · row reduction · partial pivoting · выбор ведущего элемента · число обусловленности · разреженные матрицы · fill-in · минимальная степень · законы Кирхгофа
Если каждая улица Алматы — провод, а между двумя площадями включить батарейку, ток растечётся по сотням улиц, и чтобы узнать, по каким, нужно решить систему из 7908 уравнений. Метод Гаусса «в лоб» потратил бы на неё сто шестьдесят миллиардов операций, тот же метод в удачном порядке — четверть миллиона. Здесь — почему, а ещё ведущий элемент, обусловленность и разложение, которое делают один раз.
Пусть каждая улица Алматы — провод, и сопротивление провода пропорционально длине улицы. Включим батарейку: плюс на площади Республики, минус у Зелёного базара. Ток пойдёт не одной дорогой, а растечётся по городу, и сколько его достанется каждой улице, говорят законы Кирхгофа: в каждом перекрёстке сколько втекло, столько и вытекло. Это одно линейное уравнение на перекрёсток, а перекрёстков в выгрузке OpenStreetMap, которой пользуется каталог, 7909. Ниже эта система решается в вашем браузере — методом Гаусса, но в особом порядке. Переносите полюса батарейки нажатием на карту.
Между площадью Республики и Зелёным базаром город ведёт себя как провод длиной 222 метра, хотя кратчайший путь по улицам — 3,8 км: ток делят сотни параллельных дорог. Если поставить полюс на Кок-Тобе, сопротивление подскочит до трёх с лишним километров: дорог на гору мало, и весь ток идёт по ним. А в нижней строке — цена. Метод Гаусса, который не смотрит на нули, потратил бы на эту систему 1,6·10¹¹ умножений — около минуты даже на Go, не считая полугигабайта памяти под саму матрицу. Тот же метод, исключающий перекрёстки в порядке «минимальной степени», делает четверть миллиона и укладывается в миллисекунды, а каждая новая пара полюсов стоит долю миллисекунды: тяжёлая часть работы сделана один раз.
Идея
Из первого уравнения выражаем первое неизвестное и подставляем во все остальные — неизвестных и уравнений становится на одно меньше. Повторяем, пока не останется одно уравнение с одним неизвестным; это прямой ход, и он превращает систему в лестницу. Затем обратный ход: последнее неизвестное уже известно, подставляем его в предпоследнее уравнение, и так снизу вверх. Каждое исключение — вычесть строку, умноженную на число, из остальных строк; всего около n³/3 умножений.
Подробно и с доказательствами — почему ходы не меняют решений, когда решений нет или бесконечно много и как выглядит лестница, — это разобрано в главе о методе Гаусса «Царицы наук». Здесь — инженерная сторона: откуда берутся настоящие системы, почему на диагонали оказываются нули, сколько знаков переживает ответ и как решать системы в тысячи и миллионы уравнений, где почти все коэффициенты — нули.
Уравнения цепи
Назовём φi потенциал i-го перекрёстка, а gij = 1/Rij — проводимость улицы между i и j. По закону Ома по улице течёт ток gij(φi − φj), а по закону Кирхгофа сумма вытекающих из перекрёстка токов равна току, который в него подают снаружи: +1 у плюса батарейки, −1 у минуса, 0 везде ещё. В строке i матрицы стоит сумма проводимостей всех улиц перекрёстка на диагонали и −gij напротив каждого соседа. Такую матрицу называют лапласианом графа, и у неё есть неприятность: она вырождена.
Если граф улиц связен, лапласиан L имеет ранг n − 1: решения Lφ = 0 — это ровно постоянные векторы. Если же вычеркнуть строку и столбец одного узла («заземлить» его, φ = 0), оставшаяся матрица невырождена.
Для любого вектора φ сумма φi(Lφ)i по всем узлам, если раскрыть скобки и собрать по улицам, равна Σ gij(φi − φj)² по всем улицам. Если Lφ = 0, эта сумма ноль, а все её слагаемые неотрицательны — значит, на каждой улице φi = φj, и по связности потенциал одинаков во всём городе. Постоянные векторы и правда решения: в каждой строке сумма коэффициентов ноль. Так что пространство решений одномерно, и ранг n − 1.
Пусть теперь узел z заземлён, и вектор ψ на остальных узлах решает систему без строки и столбца z. Дополним его нулём в узле z: та же сумма Σ φi(Lφ)i, взятая по узлам кроме z, по-прежнему равна сумме квадратов по улицам (слагаемое узла z нулевое, потому что φz = 0), а значит, нулю. Потенциал снова везде одинаков — и равен нулю, как в z. Единственное решение однородной системы нулевое, значит, матрица невырождена.
Физически это понятно: напряжение — разность потенциалов, и сами потенциалы определены с точностью до общей добавки, пока какой-то узел не объявлен «землёй». Если же часть цепи вообще не связана с остальной, её потенциал ничем не задан, и система вырождена даже с землёй. Симуляторы схем в таком случае сообщают о вырожденной матрице и называют «плавающий» узел.
Код: мост Уитстона
Программа решает мост Уитстона — пять резисторов и батарейку на 9 В — методом Гаусса с выбором ведущего элемента, а потом добавляет резистор, не связанный с остальной схемой. Неизвестные: ток батарейки и потенциалы трёх узлов, четвёртый узел — земля.
without row swaps: a zero pivot at step 1 with partial pivoting: potentials 9.0000 6.5138 4.1056 V battery current 47.1095 mA, the whole bridge is 191.044 ohm with a floating resistor 4-5: a zero pivot at step 6 of 6, the matrix is singular
Ведущий элемент
Первая строка вывода — не опечатка. Неизвестные записаны так: сначала ток батарейки, потом потенциалы, и первое уравнение — «φ₁ = 9 В» — тока батарейки не содержит. На диагонали ноль, делить на него нельзя, хотя система прекрасно решается. Так бывает в любой схеме с источником напряжения, если уравнения записаны в таком порядке. Лекарство — выбор ведущего элемента по столбцу: на шаге k среди оставшихся строк берут ту, где в столбце k наибольшее по модулю число, и меняют её местами с текущей.
Ноль — крайний случай. Опаснее маленькое число: формально делить можно, но множители выходят огромными, и вычитание строки, умноженной на 10¹⁴, стирает в соседних строках все знаки, кроме старших. Выбор ведущего держит множители не больше единицы по модулю, и для почти всех встречающихся на практике матриц этого хватает. Почти — потому что Джеймс Уилкинсон в 1961 году построил матрицы, на которых и с выбором ведущего числа растут вдвое на каждом шаге; в жизни такие почти не попадаются. Отдельно: для симметричных положительно определённых матриц — а лапласиан с землёй именно такой — перестановки вообще не нужны; поэтому в первой сцене строки не переставляются, а порядок выбирается только ради нулей.
Сколько знаков переживает ответ
Ошибиться можно по двум разным причинам, и их полезно не путать. Ниже решаются системы, ответ которых известен заранее — все неизвестные равны единице, — и видно, сколько верных знаков осталось.
Матрица с крошечным числом в углу — хорошая матрица: с выбором ведущего элемента ответ верен в 13–15 знаках при любом размере. Без перестановок от тех же шестнадцати знаков остаются ноль–три. Это вина алгоритма. А матрица Гильберта, с элементами 1/(i + j − 1), плоха сама по себе: её строки почти пропорциональны, и при n = 10 даже с выбором ведущего верных знаков три, при n = 12 — ни одного. Это вина задачи: малейшее изменение правой части сдвигает ответ на огромную величину, и никакой порядок действий этого не исправит. Меру этой чувствительности ввёл Алан Тьюринг в 1948 году — число обусловленности κ(A) = ‖A‖·‖A⁻¹‖.
Пусть x — точное решение Ax = b, а x̂ — приближённое, и r = b − Ax̂ — его невязка. Тогда
$$\frac{\|x - \hat x\|}{\|x\|} \le \kappa(A)\,\frac{\|r\|}{\|b\|}.$$A(x − x̂) = b − Ax̂ = r, откуда x − x̂ = A⁻¹r и ‖x − x̂‖ ≤ ‖A⁻¹‖·‖r‖. С другой стороны, ‖b‖ = ‖Ax‖ ≤ ‖A‖·‖x‖, то есть 1/‖x‖ ≤ ‖A‖/‖b‖. Перемножив два неравенства, получаем утверждение.
Метод Гаусса с выбором ведущего элемента почти всегда даёт невязку на уровне ошибок округления, около 10⁻¹⁶ от ‖b‖. Утверждение превращает её в ошибку ответа: до κ·10⁻¹⁶. Отсюда правило вычислителей: при κ ≈ 10k теряется около k знаков из шестнадцати. Маленькая невязка ещё не значит точный ответ — она значит только, что алгоритм сделал свою часть работы.
n Hilbert: kappa digits no swaps tiny corner: digits no swaps 6 2.9e+07 9 9 15 0 8 3.4e+10 6 6 14 1 10 3.5e+13 3 3 15 1 12 4.1e+16 0 0 15 0
Число обусловленности матрицы Гильберта растёт примерно в 30 раз с каждой новой строкой, и потерянные знаки идут за ним почти точно: κ = 3,4·10¹⁰ — потеряно десять знаков, остались шесть. Здесь κ посчитано по строковой норме, через точную обратную матрицу Гильберта, элементы которой — целые числа.
Разреженная матрица и заполнение
В матрице города 7909 × 7909 — 62 миллиона клеток — ненулевых всего 28 687, 0,046 %: у перекрёстка три-четыре соседа. Хранить и обрабатывать нули незачем. Но метод Гаусса их создаёт: вычитая строку k из строки i, мы пишем новые числа во все столбцы, где есть ненули в строке k. Это называют заполнением, и у него простой смысл на карте.
Исключение узла k связывает каждых двух ещё не исключённых соседей k друг с другом, и только их: новые ненули появляются ровно на месте пар соседей k, которые до этого не были связаны.
При исключении k элемент aij заменяется на aij − aikakj/akk. Поправка не равна нулю, только если aik ≠ 0 и akj ≠ 0, то есть i и j — оба соседи k. Если они уже были соседями, ненуль просто меняет значение, если нет — появляется новый. (Ещё новый ненуль мог бы случайно обратиться в ноль, но на это не рассчитывают.)
Значит, цена исключения зависит от порядка. Исключите перекрёсток в глуши, у которого два соседа, — появится одна новая связь. Исключите узел в середине, у которого к этому моменту сорок соседей, — их надо связать попарно, это сотни новых ненулей и тысячи умножений. Посмотрите, как это выглядит на куске центра 3 × 3 км.
В исходном порядке — так, как перекрёстки встретились при разборе выгрузки OpenStreetMap, — красная паутина новых связей затягивает весь район: 41 404 новых ненуля и 1,8 миллиона умножений. Обратный порядок Катхилла — Макки обходит граф слоями, в ширину, и держит ненули у диагонали узкой полосой: 7458 новых и 62 тысячи умножений. Порядок «минимальной степени» каждый раз берёт узел, у которого сейчас меньше всего соседей, и обходится 1944 новыми ненулями и 8404 умножениями. Для сравнения, плотный Гаусс на тех же 555 уравнениях — 5,7·10⁷. Найти порядок с наименьшим заполнением — задача NP-трудная (Михалис Яннакакис доказал это в 1981 году), поэтому на практике пользуются эвристиками вроде этих двух; минимальную степень описали в 1967 году Тинни и Уокер, инженеры-энергетики: у энергосистем такие же разреженные системы, только провода в них настоящие.
Разложить один раз
Ходы прямого хода выбираются по матрице и от правой части не зависят. Если записать множители, на которые вычитали строки, получится нижняя треугольная матрица L, а лестница — верхняя треугольная U, и A = LU (с перестановкой строк, если их переставляли). Решить Ax = b с готовым разложением — это два прохода по треугольным матрицам, вперёд по L и назад по U: O(n²) действий для плотной матрицы и O(число ненулей L и U) — для разреженной. Именно поэтому в первой сцене перенос полюса батарейки почти ничего не стоит: разложение не зависит от того, куда подан ток.
Программа ниже разлагает матрицу всего города один раз и решает пять систем — по одной на каждое из пяти мест, — а из этих пяти решений собирает сопротивления между всеми парами мест.
7909 intersections, 28687 nonzeros in the matrix (0.046%)
dense Gauss: about n^3/3 = 1.65e+11 multiply-adds
min degree: 36500 new nonzeros, 249163 multiply-adds
reverse Cuthill-McKee: 308386 new nonzeros, 7039994 multiply-adds
resistance between places, as a wire of this many km:
Republic Squ Green Bazaar Almaty-2 sta Botanical Ga Kok-Tobe
Republic Square 0.000 0.222 0.177 0.244 3.191
Сопротивление между двумя точками графа — это не игрушка. Оно же — «расстояние сопротивления», которым меряют, насколько хорошо связаны два места сети: чем больше независимых путей, тем оно меньше. Его же можно получить из случайных блужданий: среднее время, за которое случайный путник дойдёт из одной точки в другую и вернётся, пропорционально сопротивлению между ними.
Сколько это стоит
На шаге k прямого хода из каждой из n − k нижних строк вычитается строка длины около n − k: всего Σ(n − k)² ≈ n³/3 умножений. Обратный ход — около n²/2. Ниже замер: та же программа с тремя вложенными циклами на Python и на Go и, для сравнения, numpy.linalg.solve, которая вызывает LAPACK — то же разложение с выбором ведущего элемента, но написанное блоками под кэш процессора. Оси логарифмические, наклон прямой — показатель степени.
Наклоны близки к трём: вдвое больше уравнений — в восемь раз дольше. На сотнях уравнений Python отстаёт от LAPACK почти в тысячу раз, Go — в 7–15 раз: разница уже не в языке, а в том, как алгоритм обращается с памятью. А для системы города спасают не циклы, а нули.
| Матрица | Время | Память | Пример |
|---|---|---|---|
| плотная | ~ n³/3 | n² | тысяча уравнений — доли секунды |
| ленточная, ширина w | ~ nw² | nw | трёхдиагональная (w = 1) — O(n): прогонка |
| разреженная, хороший порядок | по заполнению | ненули L и U | 7909 перекрёстков — 2,5·10⁵ умножений |
| с готовым LU | ~ n² или ненули L и U | — | каждая новая правая часть |
Где Гаусс уже не годится — так это там, где заполнение неизбежно велико: трёхмерные сетки метода конечных элементов, где у узла десятки соседей и после любого переупорядочивания разложение не помещается в память. Там решают итерациями — методом сопряжённых градиентов и его родственниками, которые только умножают матрицу на вектор и не создают ни одного нового ненуля. Метод Гаусса там остаётся внутри — как способ построить «предобусловливатель», приближённое неполное разложение, которое ускоряет итерации.
Где это работает
- Научные библиотеки.
numpy.linalg.solveвызывает процедуру gesv из LAPACK — LU с выбором ведущего элемента; на вырожденной матрице она отвечает ошибкой «Singular matrix».scipy.sparse.linalg.spsolve— разреженное LU, по умолчанию из библиотеки SuperLU, с переупорядочиванием столбцов. - Электроника. Симуляторы схем семейства SPICE на каждом шаге по времени решают разреженную систему законов Кирхгофа — той же формы, что в первой сцене, только с транзисторами.
- Энергосистемы и инженерные расчёты. Режимы сетей, прочность конструкций, теплопроводность — разреженные системы, которые разлагают в порядке, уменьшающем заполнение.
- Внутри других методов. Неявные методы для жёстких уравнений решают систему на каждом шаге (см. стену в мороз); метод наименьших квадратов в линейной регрессии сводится к системе нормальных уравнений.
История
Самое старое известное описание метода — в восьмой главе «Девяти глав о математическом искусстве»: части этой книги написаны около 150 года до нашей эры, а в III веке Лю Хуэй снабдил её комментариями. Первая задача главы — про снопы зерна трёх сортов, а исключение ведут на счётной доске, столбцами палочек. Позже метод появляется в записях Исаака Ньютона 1669–1670 годов, где он сетует, что учебники алгебры не учат решать совместные уравнения. Гаусс в 1810 году придумал запись исключения для симметричных систем метода наименьших квадратов, и ею век с лишним пользовались профессиональные вычислители; его именем метод стали называть только в 1950-х, из-за путаницы в истории предмета.
Машинная эпоха принесла вопросы, которых у вычислителей с бумагой не было. В 1948 году Алан Тьюринг записал исключение как разложение A = LU и ввёл число обусловленности. В 1961 году Джеймс Уилкинсон показал, что метод с выбором ведущего элемента даёт точный ответ для чуть-чуть изменённой задачи, — так появился обратный анализ ошибок. А в 1960-х инженеры, считавшие энергосистемы, научились не трогать нули, и разреженные матрицы стали отдельной наукой.
Задачи
В первой задаче нужно рассчитать цепь из резисторов, в одном из тестов — кусок улиц Алматы, а в другом часть цепи висит в воздухе. Во второй метод Гаусса работает не с дробями, а в поле из двух элементов: головоломка, где каждое нажатие переключает лампочку и её соседей.
Сопротивление цепи Погасить все лампочкиДанные: Дороги Алматы © участники OpenStreetMap, ODbL