Метод Гаусса

В прошлой главе матрица превратилась в преобразование: $A$ берёт вектор $\mathbf x$ и переносит его в $A\mathbf x$. Решить систему $A\mathbf x = \mathbf b$ — значит пройти этот путь обратно и узнать, какой вектор преобразование отправило в $\mathbf b$. Для двух неизвестных хватает формулы обратной матрицы. Но у астронома, который ищет орбиту планеты, неизвестных шесть, а у инженера, который рассчитывает каркас моста, их сотни тысяч. Нужен способ, которому всё равно, сколько уравнений, и который никогда не отказывает.

Такой способ знали задолго до матриц, а знаменитым его сделала одна пропажа. В первую ночь XIX века, 1 января 1801 года, Джузеппе Пьяцци в обсерватории Палермо заметил слабую звёздочку, которой не было в каталоге. Следующей ночью она сдвинулась. Пьяцци принял её за комету, хотя для кометы она ползла слишком медленно и ровно. Он наблюдал её 24 раза, пока 11 февраля не заболел, а потом новое светило, позже названное Церерой, ушло в сияние Солнца. Чтобы осенью найти его снова, нужна была орбита, а шесть недель наблюдений давали лишь крошечный кусок пути.

Орбиту задают шесть чисел — элементы орбиты: размер и вытянутость эллипса, наклон его плоскости и линия, по которой эта плоскость пересекает плоскость земной орбиты, поворот эллипса внутри своей плоскости и положение планеты в выбранный момент. Восстановить их по наблюдениям Пьяцци взялся двадцатичетырёхлетний Карл Фридрих Гаусс. Чем кончилась охота, расскажем в конце главы. Сначала — инструмент, без которого Гаусс не обошёлся бы.

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

Три разрешённых хода

Начнём с игрушечной версии задачи астронома. Точку наблюдали трижды: в моменты $t = 1$, $2$ и $3$ она была на высоте $6$, $11$ и $18$. Допустим, высота меняется по закону $y = a + bt + ct^2$, как у брошенного камня. Каждое наблюдение даёт уравнение на неизвестные коэффициенты:

$$\left\{\begin{aligned} a + b + c &= 6, \\ a + 2b + 4c &= 11, \\ a + 3b + 9c &= 18. \end{aligned}\right.$$

Вычтем первое уравнение из второго и из третьего. Неизвестное $a$ пропадает из обоих: остаются $b + 3c = 5$ и $2b + 8c = 12$. Из последнего вычтем удвоенное предпоследнее и получим $2c = 2$. Теперь поднимаемся обратно: $c = 1$, затем $b = 5 - 3c = 2$, затем $a = 6 - b - c = 3$. Закон движения $y = 3 + 2t + t^2$, и в момент $t = 4$ точка будет на высоте $3 + 8 + 16 = 27$. Мы предсказали следующее наблюдение — в миниатюре то же, что предстояло сделать Гауссу.

Буквы в этих выкладках только занимают место: всю работу делают коэффициенты. Выпишем их таблицей, отделив чертой правые части. Каждое действие над уравнениями становится действием над строками таблицы:

$$\left(\begin{array}{ccc|c} 1 & 1 & 1 & 6 \\ 1 & 2 & 4 & 11 \\ 1 & 3 & 9 & 18 \end{array}\right) \xrightarrow{\substack{R_2 - R_1 \\ R_3 - R_1}} \left(\begin{array}{ccc|c} 1 & 1 & 1 & 6 \\ 0 & 1 & 3 & 5 \\ 0 & 2 & 8 & 12 \end{array}\right) \xrightarrow{R_3 - 2R_2} \left(\begin{array}{ccc|c} 1 & 1 & 1 & 6 \\ 0 & 1 & 3 & 5 \\ 0 & 0 & 2 & 2 \end{array}\right).$$

Над уравнениями ходы делают то же, что мы делали руками: переставляют уравнения, умножают уравнение на число, прибавляют к одному уравнению другое. Но почему после таких ходов решения не теряются и не появляются лишние? Это нужно доказать, и доказательство здесь — самое наглядное место всего метода.

Хитрость в том, что каждый ход можно отменить другим ходом. Поэтому достаточно проверить одно: после хода старые решения остаются решениями.

Начнём с двух уравнений с двумя неизвестными. На чертеже это две прямые, $R_1$ и $R_2$ (сначала $x + y = 3$ и $x - y = 1$), а их общая точка $P$ — решение системы. Заменим второе уравнение на $R_2 + kR_1$. Точка $P$ обращает оба уравнения в верные равенства. Сложим эти равенства, умножив первое на $k$, — получим верное равенство, то есть $P$ удовлетворяет и новому уравнению. Новая прямая при любом $k$ проходит через $P$: меняя $k$, мы поворачиваем вторую прямую вокруг решения. Такое семейство прямых через одну точку называют пучком. Подберём $k$ так, чтобы из уравнения пропал $x$ (сначала это $k = -1$: $(x - y) - (x + y) = 1 - 3$, то есть $-2y = -2$). Прямая стала горизонтальной и по-прежнему проходит через $P$. Вот что такое исключение неизвестного на чертеже: поворот прямой вокруг решения, пока она не ляжет параллельно оси. В системе из любого числа уравнений рассуждение то же. Ход $R_i \to R_i + kR_j$ меняет только уравнение $i$, и если $\mathbf x$ удовлетворял уравнениям $i$ и $j$, то удовлетворяет и новому: левая часть нового уравнения — сумма левых частей с множителями $1$ и $k$, правая — сумма правых с теми же множителями. Умножение уравнения на число и перестановка уравнений тем более сохраняют решения. Значит, каждое решение старой системы — решение новой. Обратно: ход $R_i \to R_i + kR_j$ отменяется ходом $R_i \to R_i - kR_j$ (строка $R_j$ осталась прежней), умножение на $k \ne 0$ — умножением на $\frac1k$, перестановка — той же перестановкой. Применим к новой системе отменяющий ход: по уже доказанному её решения останутся решениями старой. У обеих систем одни и те же решения. Доведём исключение до конца, и прямые станут вертикальной и горизонтальной, $x = 2$ и $y = 1$: ответ читается сразу. Заодно понятно, почему нельзя умножать строку на ноль. Такой ход не отменить, а получившемуся уравнению $0 = 0$ удовлетворяет любая точка.

В системе два уравнения: $x + y = 3$ и $x - y = 1$. Второе уравнение умножили на $0$. Что стало с множеством решений?

Не изменилось: умножение — разрешённый ход Стало бесконечным: осталась одна прямая $x + y = 3$ Стало пустым: $0 = 1$ Пустым, потому что пропали оба неизвестных

Уравнение $0x + 0y = 0$ верно при любых $x$ и $y$, так что от системы осталось одно уравнение $x + y = 3$ — целая прямая решений. Умножать можно только на ненулевое число: такой ход обратим, и решения не меняются.

Лестница

Ходы есть — нужна стратегия. Цель — матрица, в которой под «лестницей» стоят одни нули, а каждая следующая ступенька начинается правее предыдущей:

$$\left(\begin{array}{cccc|c} \boxed{2} & 1 & 0 & 3 & 5 \\ 0 & 0 & \boxed{-1} & 4 & 1 \\ 0 & 0 & 0 & \boxed{7} & 14 \\ 0 & 0 & 0 & 0 & 0 \end{array}\right)$$

Ступенчатую систему решают снизу вверх, как мы уже делали: последнее ненулевое уравнение содержит меньше всего неизвестных, его решаем первым и подставляем найденное выше. Это обратный ход. А получают лестницу прямым ходом. Идём по столбцам слева направо. В очередном столбце ищем ненулевой элемент среди строк, которые ещё не стали ступеньками; если его нет, столбец пропускаем. Найденный элемент перестановкой поднимаем на верх оставшейся части — это новый ведущий. Из каждой строки ниже вычитаем ведущую строку с таким множителем, чтобы под ведущим получился ноль: если ведущий равен $p$, а под ним стоит $q$, множитель равен $\frac qp$. Затем опускаемся на строку ниже и переходим к следующему столбцу.

Доказываем индукцией по числу столбцов: прямой ход устраивает первый столбец и дальше работает с матрицей поменьше.

Если столбец один, ход делает так: если в столбце есть ненулевой элемент, поднимает его наверх и обнуляет всё под ним, а если все элементы нули, ничего не делает. В обоих случаях получается ступенчатый столбец. Пусть для матриц с меньшим числом столбцов утверждение доказано. Если первый столбец целиком нулевой, отбросим его: к остальной матрице ход применяется так же, и по предположению индукции она станет ступенчатой. Нулевой столбец слева лестницу не портит — все ступеньки просто сдвинуты на столбец вправо. Если в первом столбце есть ненулевой элемент $p$, перестановка поднимает его в первую строку, а вычитания делают нулями все элементы под ним: из строки с элементом $q$ вычитаем первую строку, умноженную на $\frac qp$, и на месте $q$ остаётся $q - \frac qp \cdot p = 0$. Дальше ход не трогает ни первую строку, ни первый столбец (там уже нули под $p$), а работает с матрицей из остальных строк без первого столбца. По предположению индукции она станет ступенчатой. Сверху стоит строка с ведущим $p$ в первом столбце, под ней — ступенчатая матрица, у которой первый столбец нулевой. Ступеньки нижней части начинаются правее первого столбца, так что вся матрица ступенчатая.

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

$$\left(\begin{array}{ccc|c} 1 & 1 & 1 & 6 \\ 0 & 1 & 3 & 5 \\ 0 & 0 & 2 & 2 \end{array}\right) \longrightarrow \left(\begin{array}{ccc|c} 1 & 0 & 0 & 3 \\ 0 & 1 & 0 & 2 \\ 0 & 0 & 1 & 1 \end{array}\right),$$

и ответ $a = 3$, $b = 2$, $c = 1$ стоит прямо в последнем столбце.

Назван он по Вильгельму Йордану, немецкому геодезисту, который описал этот вариант в учебнике геодезии в конце XIX века. По-русски его фамилию обычно пишут на французский лад, и метод нередко по ошибке связывают с Камилем Жорданом, чья нормальная форма матриц встретится в главе 38.

Иногда на месте очередного ведущего оказывается ноль. Возьмём систему $x + y + 2z = 5$, $2x + 2y + 5z = 11$, $3x + 4y + 7z = 17$. После вычитаний $R_2 - 2R_1$ и $R_3 - 3R_1$ из второго уравнения пропали сразу и $x$, и $y$: оно стало $z = 1$. Ведущий второго столбца надо искать ниже — он в третьей строке, $y + z = 2$. Меняем строки местами, и лестница готова. Попробуйте провести это сами.

Решите систему $x + y + z = 4$, $2x + 3y + z = 4$, $x - y + 2z = 9$. Ответ запишите тройкой $(x;\,y;\,z)$.

$R_2 - 2R_1$ даёт $y - z = -4$, $R_3 - R_1$ даёт $-2y + z = 5$. Прибавим к третьему удвоенное второе: $-z = -3$, то есть $z = 3$. Тогда $y = -4 + z = -1$, а $x = 4 - y - z = 2$. Проверка по третьему уравнению: $2 + 1 + 6 = 9$.

Одно, бесконечно много или ни одного

До сих пор лестница выходила полной: сколько неизвестных, столько и ступенек. Так бывает не всегда. Возьмём систему

$$\left\{\begin{aligned} x + 2y + z &= 4, \\ 2x + 4y + 3z &= 9, \\ 3x + 6y + 4z &= 13. \end{aligned}\right.$$

Вычитания $R_2 - 2R_1$ и $R_3 - 3R_1$ дают одинаковые строки $z = 1$, и после $R_3 - R_2$ третье уравнение превращается в $0 = 0$. Обратный ход приводит к приведённому ступенчатому виду:

$$\left(\begin{array}{ccc|c} 1 & 2 & 0 & 3 \\ 0 & 0 & 1 & 1 \\ 0 & 0 & 0 & 0 \end{array}\right), \qquad\text{то есть}\qquad x + 2y = 3, \quad z = 1.$$

Ступенька во втором столбце так и не появилась, и неизвестное $y$ ничем не ограничено. Задайте его как угодно — остальные подстроятся: $x = 3 - 2y$, $z = 1$. Решений бесконечно много, и все они лежат на одной прямой в пространстве.

Любое решение системы — столбец $(x;\,y;\,z)$. Частное решение: свободному неизвестному дали значение $0$, главные нашли обратным ходом. Параметр — значение свободного неизвестного $y$. Он пробегает все действительные числа. Направление прямой решений. Этот столбец решает однородную систему с нулевыми правыми частями: сдвиг вдоль него не меняет левых частей уравнений. Пример: при $t = 2$ получаем $(-1;\,2;\,1)$. Проверка: $-1 + 4 + 1 = 4$, $-2 + 8 + 3 = 9$, $-3 + 12 + 4 = 13$.

Если свободных неизвестных два, в общем решении два параметра, и решения заполняют плоскость. А если в третьем уравнении справа стоит не $13$, а $14$, та же цепочка ходов закончится строкой $0 = 1$. Это уравнение неверно при любых $x$, $y$, $z$, и система несовместна.

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

Идея: приведённый вид можно прочитать по связям между столбцами матрицы, а эти связи ходы не меняют.

Обозначим столбцы матрицы $A$ через $\mathbf a_1, \dots, \mathbf a_n$. Равенство $c_1\mathbf a_1 + \dots + c_n\mathbf a_n = \mathbf 0$ означает ровно то, что столбец $(c_1;\,\dots;\,c_n)$ решает однородную систему $A\mathbf x = \mathbf 0$: $i$-я строка этого векторного равенства и есть $i$-е уравнение системы. Если матрица $B$ получена из $A$ элементарными преобразованиями строк, то по предыдущей теореме у систем $A\mathbf x = \mathbf 0$ и $B\mathbf x = \mathbf 0$ одни и те же решения. Значит, любая линейная комбинация столбцов $A$ равна нулю тогда и только тогда, когда та же комбинация столбцов $B$ равна нулю. В частности, столбец $A$ выражается через предыдущие столбцы с некоторыми коэффициентами ровно тогда, когда тот же столбец $B$ выражается через предыдущие столбцы $B$ с теми же коэффициентами. Посмотрим на приведённый ступенчатый вид $R$. Столбец с $k$-м по счёту ведущим элементом — это $k$-й единичный столбец $\mathbf e_k$: единица в строке $k$, остальные нули. Он не выражается через предыдущие столбцы, потому что у всех них в строке $k$ стоят нули. Столбец без ведущего имеет ненулевые элементы только в строках уже встреченных ведущих, поэтому он равен $r_1\mathbf e_1 + \dots + r_k\mathbf e_k$, где $r_1, \dots, r_k$ — его элементы, и выражается через предыдущие ведущие столбцы. Выходит, ведущие столбцы $R$ — ровно те, что не выражаются через предыдущие. По второму шагу это свойство у столбца $R$ такое же, как у столбца $A$ с тем же номером, так что положения ступенек определены самой матрицей $A$. Столбец $R$ без ведущего равен $r_1\mathbf e_1 + \dots + r_k\mathbf e_k$, и такая запись единственна: коэффициент при $\mathbf e_i$ — просто $i$-й элемент столбца. По второму шагу соответствующий столбец $A$ выражается через ведущие столбцы $A$ с теми же коэффициентами $r_1, \dots, r_k$. Если бы через ведущие столбцы $A$ он выражался двумя способами, то и через $\mathbf e_1, \dots, \mathbf e_k$ — тоже двумя, а это невозможно. Значит, числа $r_i$ определены матрицей $A$. Положения ведущих единиц и все остальные элементы приведённого вида определены матрицей $A$, а не выбором ходов. Два приведённых вида одной матрицы совпадают.

Любой ступенчатый вид дополировывается до приведённого, не меняя числа ступенек, поэтому ранг можно считать по любой лестнице. У системы про параболу $\operatorname{rank} A = 3$, у системы со свободным $y$ ранг равен $2$. Теперь можно сказать, когда у системы есть решения, одним сравнением двух чисел.

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

Приведём $(A \mid \mathbf b)$ к ступенчатому виду. Без последнего столбца получится ступенчатый вид самой $A$: все ходы выбирались по столбцам слева направо. Поэтому ранги различаются ровно тогда, когда ведущий элемент оказался в последнем столбце, то есть появилась строка $(0\ \dots\ 0 \mid c)$ с $c \ne 0$. На чертеже три уравнения с двумя неизвестными — три прямые. Первые две пересекаются в точке $P$ и дают ступеньки в столбцах $x$ и $y$. Прямой ход вычитает из третьего уравнения комбинацию $\alpha R_1 + \beta R_2$ с той же левой частью. Все комбинации $R_1$ и $R_2$ — прямые пучка через $P$, а нужная из них — прямая через $P$, параллельная третьей. Если третья прямая проходит мимо $P$, от третьего уравнения остаётся $0 \cdot x + 0 \cdot y = c$ с $c \ne 0$. Этому уравнению не удовлетворяет ни одна точка: решений нет, и $\operatorname{rank}\,(A \mid \mathbf b) = 3 > 2 = \operatorname{rank} A$. Чем дальше третья прямая от $P$, тем больше $|c|$. Если ведущего в последнем столбце нет, все ступеньки стоят в столбцах неизвестных. Дадим свободным неизвестным любые значения и найдём главные обратным ходом: каждое ненулевое уравнение выражает своё главное неизвестное через те, что правее. Решение существует. Протащите третью прямую через $P$ — строка обнулится целиком, $0 = 0$, и $P$ станет решением всей системы. Система совместна ровно тогда, когда ранги равны. Если при этом $r = n$, свободных неизвестных нет, обратный ход определяет каждое неизвестное однозначно, и решение одно. Если $r < n$, свободных неизвестных $n - r$; каждое принимает любое значение, и разные значения дают разные решения — их бесконечно много.

В разных странах эту теорему называют по-разному. У нас и в Польше — по Леопольду Кронекеру и Альфредо Капелли, в Италии и в англоязычных книгах — теоремой Руше — Капелли, во Франции — Руше — Фонтене, в Испании — Руше — Фробениуса.

В пространстве каждое уравнение с тремя неизвестными задаёт плоскость, и теорема превращается в каталог того, как могут расположиться три плоскости:

Строка «книжка» включает и случай, когда две плоскости совпали, а третья их пересекает. «Призма» — это три плоскости, которые попарно пересекаются по трём параллельным прямым; другой несовместный случай того же ранга — две параллельные плоскости и третья, секущая обе. В последней строке параллельны все три, и хотя бы две из них различны.

Особенно простой случай — однородная система, у которой все правые части нули. Она всегда совместна: нулевой столбец её решает. Интересно, есть ли у неё другие решения.

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

В системе три уравнения и четыре неизвестных. Может ли у неё быть ровно одно решение?

Да, если уравнения независимы Нет: решений либо нет, либо бесконечно много Нет: решений всегда бесконечно много

Ранг не больше числа строк, то есть $3 < 4$, так что свободное неизвестное есть всегда. Если система совместна, оно принимает любые значения. Но совместной она может и не быть: например, если два уравнения отличаются только правой частью.

Найдите ранг матрицы $\begin{pmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{pmatrix}$.

$R_2 - 4R_1$ и $R_3 - 7R_1$ дают строки $(0;\,-3;\,-6)$ и $(0;\,-6;\,-12)$. Вторая из них — удвоенная первая, и после $R_3 - 2R_2$ остаётся нулевая строка. Ступенек две: ранг равен $2$.

При каком $a$ система $x + y + z = 1$, $x + 2y + 3z = 2$, $2x + 3y + (a + 6)z = 4$ не имеет решений?

$R_2 - R_1$ даёт $y + 2z = 1$, $R_3 - 2R_1$ даёт $y + (a + 4)z = 2$. Вычтем одно из другого: $(a + 2)z = 1$. При $a \ne -2$ отсюда находится $z$, и решение единственно. При $a = -2$ получается $0 = 1$: ранг матрицы $2$, ранг расширенной $3$, решений нет.

Сколько стоит решение

Гаусс считал вручную, и каждое умножение отнимало время. Сколько действий нужно, чтобы решить методом Гаусса систему из $n$ уравнений с $n$ неизвестными?

Хитрость в том, чтобы увидеть, где работа сосредоточена: в квадратах, которые с каждым шагом уменьшаются.

Первый шаг прямого хода. Под ведущим $n - 1$ строк. Для каждой нужен множитель (одно деление) и обновление остальных $n$ чисел строки: $n - 1$ элемент матрицы и правая часть. Каждое обновление — одно умножение и одно вычитание. Основная работа идёт в квадрате $(n - 1) \times (n - 1)$ в правом нижнем углу матрицы. Второй шаг работает с квадратом $(n - 2) \times (n - 2)$, третий — ещё меньше. Квадраты вложены друг в друга и сжимаются к углу. Всего в них $(n - 1)^2 + (n - 2)^2 + \dots + 1^2$ клеток, и каждая клетка — умножение плюс вычитание. Поставим эти числа столбиками ширины $1$: столбик высоты $m^2$ над отрезком $[m;\,m + 1]$ целиком лежит под параболой $y = x^2$, а тот же столбик над отрезком $[m - 1;\,m]$ целиком её накрывает. Значит, сумма столбиков зажата между площадями под параболой: $\frac{(n - 1)^3}{3} \le 1^2 + 2^2 + \dots + (n - 1)^2 \le \frac{n^3}{3}$. Площадь под параболой мы считали в главе об интеграле. Точный счёт. На шаге, где под ведущим $m$ строк, делается $m(m + 2)$ умножений и делений и $m(m + 1)$ вычитаний. Обратный ход: неизвестное номер $i$ стоит $n - i$ умножений, столько же вычитаний и одно деление. Складывая с формулами $1 + 2 + \dots + k = \frac{k(k + 1)}{2}$ и $1^2 + \dots + k^2 = \frac{k(k + 1)(2k + 1)}{6}$ из главы о последовательностях, получаем $\sum_{m=1}^{n-1} m(m + 2) + \frac{n(n + 1)}{2} = \frac{n^3}{3} + n^2 - \frac n3$ и $\sum_{m=1}^{n-1} m(m + 1) + \frac{n(n - 1)}{2} = \frac{n^3}{3} + \frac{n^2}{2} - \frac{5n}{6}$. Всё, кроме квадратов, — множители, правые части, обратный ход — добавляет лишь слагаемые порядка $n^2$. При больших $n$ главное — $\frac{2n^3}{3}$ действий: в два раза больше неизвестных — в восемь раз больше работы. Главная часть: обновление сжимающихся квадратов при прямом ходе, по умножению и вычитанию на клетку. Множители, правые части и обратный ход — по нескольку действий на каждый элемент матрицы. Поправка, которая остаётся от точных сумм; при больших $n$ она ничтожна. Пример: при $n = 3$ получаем $18 + 13{,}5 - 3{,}5 = 28$ действий ($17$ умножений и делений и $11$ сложений и вычитаний). Для шести неизвестных Гаусса — $191$ действие, для тысячи — $668\,165\,500$, и из них $666\,666\,667$ даёт первое слагаемое.

Много это или мало? Сравните с формулами Крамера из главы о системах: каждое неизвестное — отношение двух определителей. Если считать определитель раскрытием по строке, в нём $n!$ слагаемых. Для $n = 20$ это около $10^{21}$ умножений: при миллиарде операций в секунду — тридцать тысяч лет. Метод Гаусса решает ту же систему за $5910$ действий, то есть за несколько микросекунд.

Кубический рост — главный закон этого раздела. Удвоили размер сетки в инженерной модели — работа выросла в восемь раз. Поэтому для огромных систем, где почти все коэффициенты нули, придумано множество обходных путей, которые пользуются этими нулями. Но плотную систему без особой структуры по-прежнему решают методом Гаусса. Им же измеряют мощь суперкомпьютеров: в рейтинге TOP500, который ведут с 1993 года, машины соревнуются в одном упражнении — решить огромную плотную систему линейных уравнений.

Сколько умножений и делений понадобится методу Гаусса для системы из $10$ уравнений с $10$ неизвестными (с единственным решением и без перестановок)?

По формуле $\frac{n^3}{3} + n^2 - \frac n3$ при $n = 10$: $\frac{1000}{3} + 100 - \frac{10}{3} = \frac{990}{3} + 100 = 330 + 100 = 430$.

Тонкий лёд

Компьютер хранит числа с конечным числом значащих цифр — обычно около шестнадцати. После каждого действия результат округляется, и ошибки накапливаются. Для метода Гаусса это не мелочь: неудачный ведущий элемент может погубить весь ответ. Покажем это на калькуляторе, который помнит только три значащие цифры. Решим систему

$$\left\{\begin{aligned} 0{,}0001x + y &= 1, \\ x + y &= 2. \end{aligned}\right.$$

Её точное решение $x = \frac{10\,000}{9999} \approx 1{,}0001$, $y = \frac{9998}{9999} \approx 0{,}9999$. Ведущий — $0{,}0001$, множитель $10\,000$. Вычитаем: $(1 - 10\,000)y = 2 - 10\,000$, то есть $-9999y = -9998$. Три цифры калькулятора превращают обе части в $-10\,000$, и выходит $y = 1$. Подставляем в первое уравнение: $0{,}0001x = 1 - 1 = 0$, то есть $x = 0$. Ответ неверен целиком. Двойка из второго уравнения утонула в огромном $10\,000$, а потом крошечный ведущий раздул оставшуюся ошибку.

Поменяем уравнения местами. Ведущий теперь $1$, множитель $0{,}0001$: $(1 - 0{,}0001)y = 1 - 0{,}0002$, после округления $y = 1$, и $x = 2 - y = 1$. Ответ верен с точностью до третьей цифры.

Выбор главного элемента исправляет ошибки алгоритма, но не может исправить саму задачу. Вспомните систему из главы о системах: $x + y = 2$, $x + 1{,}01y = 2{,}01$ с решением $(1;\,1)$. Сдвиньте правую часть на сотую, до $2{,}02$, и решение уедет в $(0;\,2)$. Правая часть изменилась на полпроцента, а решение — на сто процентов. Прямые почти параллельны, и никакой порядок действий этого не исправит: так устроены сами данные.

У матрицы из нашего примера $\kappa \approx 402$, и полпроцента превратились в сотню процентов — в двести раз больше, в пределах обещанного. Грубое правило вычислителей: при $\kappa \approx 10^k$ в ответе теряется около $k$ верных знаков. У матрицы Гильберта $10 \times 10$ с элементами $\frac{1}{i + j - 1}$ число обусловленности около $1{,}6 \cdot 10^{13}$, и из шестнадцати знаков компьютера верными останутся лишь два-три.

Записать ходы: LU-разложение

Астроном получает новые наблюдения, а моменты времени, в которые он их делает, те же. В нашей игрушке это значит: матрица $A$ та же, меняются только правые части. Ходы прямого хода выбираются по матрице и от правых частей не зависят. Их стоит записать один раз.

$L$ (от lower — «нижняя»): журнал ходов. Строки $2$ и $3$ получили $R_1$ с множителем $1$, а потом из строки $3$ вычли $R_2$ с множителем $2$. $U$ (от upper — «верхняя»): лестница, которой закончился прямой ход в задаче про параболу. Пример: новые наблюдения $2$, $4$, $8$ в те же моменты. Решаем $L\mathbf z = \mathbf b$ сверху вниз: $z_1 = 2$, $z_2 = 4 - 2 = 2$, $z_3 = 8 - 2 - 2 \cdot 2 = 2$. Затем $U\mathbf x = \mathbf z$ снизу вверх: $c = 1$, $b = 2 - 3 = -1$, $a = 2 + 1 - 1 = 2$. Закон $y = 2 - t + t^2$ проходит через все три точки.

Две треугольные системы решаются примерно за $2n^2$ действий против $\frac23 n^3$ у полного метода. Для тысячи неизвестных это два миллиона действий вместо почти семисот миллионов. Если перестановки нужны, разложение записывают как $PA = LU$: матрица $P$ переставляет строки $A$ так, как это сделал бы выбор главного элемента, после чего прямой ход идёт без перестановок. Именно в таком виде метод Гаусса живёт в вычислительных библиотеках. LU-разложение описал в 1938 году польский астроном Тадеуш Банахевич. В 1948 году Алан Тьюринг разобрал, как в нём накапливаются ошибки округления, и ввёл термин «число обусловленности».

Церера найдена

Вернёмся в 1801 год. Гаусс придумал, как восстановить орбиту по трём наблюдениям, вычислил, где искать Цереру, и отправил результат астроному Францу Ксаверу фон Цаху. 31 декабря 1801 года фон Цах, а следующей ночью и Генрих Ольберс нашли её меньше чем в полуградусе от места, указанного Гауссом. Двадцатичетырёхлетний математик в одночасье стал знаменитым астрономом. Через несколько лет он возглавил обсерваторию в Гёттингене.

Свой способ он изложил в 1809 году в книге «Теория движения небесных тел, обращающихся вокруг Солнца по коническим сечениям». Там же — метод наименьших квадратов: как согласовать десятки наблюдений, каждое со своей погрешностью. О нём и о споре Гаусса с Лежандром за первенство расскажет глава 39. Метод наименьших квадратов приводит к системам линейных уравнений, и Гаусс объяснил, как решать их последовательным исключением. В 1810 году, рассчитывая орбиту Паллады — второго астероида, открытого Ольберсом весной 1802-го, — он придумал для исключения компактную запись. Профессиональные вычислители пользовались ею весь XIX век.

Сам способ старше Гаусса почти на два тысячелетия: исключение неизвестных описано в восьмой главе китайской «Математики в девяти книгах», и счётную доску со столбцами палочек мы видели в главе о системах. В Европе ему учили по «Всеобщей арифметике» Исаака Ньютона, изданной в 1707 году по его лекциям. К концу XVIII века исключение неизвестных стало обычным уроком алгебры. Имя Гаусса прикрепилось к школьному алгоритму лишь в 1950-х годах, в эпоху первых компьютеров, — во многом по путанице в его истории.

Сегодня Церера — карликовая планета, самое крупное тело пояса астероидов между Марсом и Юпитером. В 2015–2018 годах вокруг неё летал зонд Dawn и фотографировал её яркие пятна — отложения солей. Траекторию зонда уточняли по измерениям методом наименьших квадратов, прямым потомком расчётов Гаусса.

Куда дальше

В этой главе складывались очень разные вещи. Строки матрицы мы прибавляли одну к другой, как векторы. Решения однородной системы можно сложить или умножить на число — и снова получится решение. Параболы $y = a + bt + ct^2$ складываются коэффициентами: сумма парабол из наших двух примеров, $3 + 2t + t^2$ и $2 - t + t^2$, снова многочлен той же формы, и он проходит через точки, высоты которых равны суммам высот. Звуковые сигналы складываются, когда играют два инструмента, и умножаются на число, когда крутят ручку громкости. Стрелки, строки, многочлены, сигналы — почему все они ведут себя одинаково, и что можно понять сразу про все такие миры? Ответ — одно определение на восемь строк. С него начинается глава 37.