Матричный метод решения систем ⚖️

У тебя уже есть два рабочих инструмента для системы линейных уравнений: метод Гаусса, который перемалывает любую систему в ступенчатый вид, и правило Крамера, которое выдаёт каждую неизвестную отдельной дробью из определителей. Есть и третий способ, и он самый короткий из всех — настолько короткий, что записывается в одну строку:

Всё. Система из ста уравнений со ста неизвестными сворачивается в одну формулу, где нет ни ступенчатого вида, ни ведущих элементов, ни ста первых определителей. Выглядит как победа. Ровно поэтому матричный метод любят учебники: он красиво замыкает предыдущую тему про обратную матрицу и превращает решение системы в обычную алгебру, где $A$ можно «перенести в другую часть», поделив на неё.

А теперь неприятная правда, которую редко проговаривают вслух: в реальном численном коде так почти никогда не пишут. Если ты откроешь исходники любой серьёзной библиотеки — NumPy, SciPy, PyTorch, Eigen, LAPACK, — ты не найдёшь там строчки «посчитать обратную матрицу и умножить на правую часть» в качестве способа решить систему. Более того, в кодстайлах ML-команд, в книгах по численным методам и в документации самой NumPy прямым текстом написано: не вычисляй inv(A), если тебе нужен только $x$. Это не суеверие и не вкусовщина — за этим стоят конкретные цифры по точности, скорости и памяти, которые мы в этом уроке измерим своими руками.

Так что этот урок — про две вещи сразу. Первая: как матричный метод устроен, когда он законен, как им пользоваться и в каком редком, но реальном сценарии он действительно окупается (спойлер: когда правых частей много, а матрица одна). Вторая — и она важнее — про то, насколько вообще можно доверять ответу, который выдаёт решатель систем. Здесь появится центральное понятие всей вычислительной линейной алгебры: число обусловленности матрицы. Это одно число, которое отвечает на вопрос «если мои входные данные известны с точностью до четвёртого знака, со сколькими верными знаками я получу ответ?». Оказывается, бывают абсолютно нормальные с виду системы $2\times2$ из целых чисел, где изменение правой части на одну сотую процента меняет ответ на 400 процентов. И ты увидишь эти числа посчитанными, а не рассказанными.

Понимание обусловленности — это, пожалуй, самый практичный навык из всей линейной алгебры для того, кто идёт в ML. Именно оно объясняет, почему коэффициенты регрессии скачут при добавлении одного наблюдения, почему нейросеть не сходится без нормализации признаков, почему в гауссовских процессах к ковариационной матрице зачем-то прибавляют крошечное число на диагональ, и почему две математически эквивалентные строчки кода дают ответы, различающиеся в восьмом знаке — или в первом.

🎯 Ты узнаешь:

• как вывести формулу $x = A^{-1}b$ и почему домножать нужно строго слева, а справа — нельзя даже формально;

• как решить систему $3\times3$ матричным методом от начала до конца и проверить ответ подстановкой;

• что такое система с несколькими правыми частями $AX = B$ и когда явная обратная матрица действительно окупается;

• честное сравнение трёх методов — Гаусса, Крамера и матричного — по числу операций, устойчивости и области применимости, с конкретными цифрами;

• что такое норма вектора и норма матрицы, как считать их руками и что они измеряют;

• что такое число обусловленности $\operatorname{cond}(A) = \|A\| \cdot \|A^{-1}\|$ и как по нему предсказать, сколько верных знаков останется в ответе;

• почему в коде пишут numpy.linalg.solve(A, b), а не numpy.linalg.inv(A) @ b — с измеренными числами по точности и по времени;

• чем «почти вырожденная» матрица опаснее честно вырожденной: одна падает с ошибкой, вторая тихо возвращает мусор.

История: откуда это взялось?

Сама идея «перенести матрицу в другую часть уравнения» появилась вместе с понятием обратной матрицы — у Артура Кэли в мемуаре 1858 года. Кэли показал, что матрицы образуют алгебру, в которой можно складывать, умножать и (иногда) обращать, и запись $x = A^{-1}b$ была для него естественным следствием: если $A$ обратима, то уравнение $Ax = b$ решается ровно так же, как школьное $ax = b$ решается делением на $a$. Для XIX века это было изящное концептуальное завершение теории.

Но уже к середине XX века, когда появились первые электронные вычислительные машины и на них стали считать реальные инженерные задачи, обнаружилось нечто неожиданное: формула, безупречная на бумаге, на машине ведёт себя отвратительно. В 1947 году Джон фон Нейман и Герман Голдстайн опубликовали работу «Numerical Inverting of Matrices of High Order» — первое серьёзное исследование того, что происходит с точностью, когда матрицу обращают в арифметике с конечным числом разрядов. Их вывод был отрезвляющим: ошибка округления не просто накапливается, а усиливается, и коэффициент усиления зависит от свойств самой матрицы.

Через год, в 1948-м, Алан Тьюринг — тот самый Тьюринг — опубликовал статью «Rounding-off Errors in Matrix Processes», где ввёл в оборот сам термин condition number, число обусловленности. Тьюринг занимался тогда проектированием вычислительной машины ACE и на практике столкнулся с тем, что одни системы решаются надёжно, а другие разваливаются. Он предложил измерять «трудность» матрицы одним числом и показал, что именно оно управляет тем, насколько ошибка входных данных раздувается в ответе. Занятно, что понятие, без которого сегодня не обходится ни один курс численных методов, родилось не в чистой математике, а из инженерной необходимости отладить конкретный компьютер.

Довёл теорию до законченного вида англичанин Джеймс Уилкинсон в 1950–1960-х годах. Он разработал обратный анализ ошибок (backward error analysis): вместо того чтобы спрашивать «насколько неточен полученный ответ», Уилкинсон предложил спрашивать «какую задачу мы на самом деле решили точно». Оказалось, что метод Гаусса с выбором ведущего элемента решает точно систему с чуть-чуть возмущённой матрицей — и весь вопрос сводится к тому, насколько сильно решение реагирует на это «чуть-чуть». А реакция как раз и измеряется числом обусловленности. За эти работы Уилкинсон получил премию Тьюринга в 1970 году, а его книга «The Algebraic Eigenvalue Problem» (1965) до сих пор считается классикой.

Отдельный сюжет — матрица Гильберта. В 1894 году Давид Гильберт, занимаясь задачей приближения функций многочленами, наткнулся на матрицу с элементами $h_{ij} = \dfrac{1}{i+j-1}$. Она симметричная, все элементы положительные и вполне «обычные» на вид, определитель не равен нулю — то есть формально всё в порядке. Но вычислительно эта матрица оказалась чудовищем: уже для размера $10\times10$ её число обусловленности превышает $10^{13}$, и решить с ней систему в обычной машинной точности практически невозможно. Гильберт нашёл её задолго до появления компьютеров, но с 1950-х она стала стандартным тестовым полигоном для численных алгоритмов — и мы тоже прогоним её через все наши формулы в этом уроке.

Практический итог этой истории такой. Библиотека LINPACK (1979), а затем LAPACK (1992), на которых сегодня стоит вся вычислительная линейная алгебра, включая NumPy и SciPy, устроены так, что функции «решить систему» и «обратить матрицу» — это разные функции с разными алгоритмами, и первая никогда не реализуется через вторую. Это прямое наследие фон Неймана, Тьюринга и Уилкинсона: они выяснили, что решать систему через явную обратную матрицу — это делать лишнюю работу и одновременно портить точность.

Матричный метод: вывод формулы

Интуиция: система как одно уравнение

Возьми обычную систему трёх уравнений с тремя неизвестными:

Ты уже умеешь записывать её в матричном виде $Ax = b$, где

Посмотри на запись $Ax = b$ и на школьное уравнение $5x = 20$. Они выглядят одинаково. В школьном случае ты делишь обе части на 5 и получаешь $x = 20/5 = 4$. Соблазн сделать то же самое с матрицами огромен — и в этом соблазне есть здоровое зерно. Деления на матрицу не существует, но существует умножение на обратную матрицу, а это ровно то же самое по смыслу: $A^{-1}$ отменяет действие $A$, как $1/5$ отменяет умножение на 5.

Разница только одна, зато принципиальная: умножение матриц некоммутативно. У числа $1/5$ нет «стороны» — умножать на него можно хоть слева, хоть справа, результат один. У матрицы $A^{-1}$ сторона есть, и выбор стороны определяет, получится у тебя решение или бессмыслица.

Вывод

Пусть дана система $Ax = b$, где $A$ — квадратная матрица порядка $n$, $x$ и $b$ — векторы-столбцы из $n$ элементов. Предположим, что $\det A \ne 0$, то есть $A$ обратима.

Домножим обе части равенства на $A^{-1}$ слева:

Слева воспользуемся ассоциативностью матричного умножения (переставлять множители нельзя, а группировать — можно):

По определению обратной матрицы $A^{-1}A = E$, а умножение на единичную матрицу ничего не меняет:

Матричный метод (формула решения): если $A$ — квадратная матрица порядка $n$ и $\det A \ne 0$, то система $Ax = b$ имеет решение

$$x = A^{-1}b$$

и это решение единственно.

Единственность здесь получается бесплатно, и это важно проговорить. Допустим, у системы есть два решения $x_1$ и $x_2$, то есть $Ax_1 = b$ и $Ax_2 = b$. Домножим оба равенства слева на $A^{-1}$: получим $x_1 = A^{-1}b$ и $x_2 = A^{-1}b$, то есть $x_1 = x_2$. Никаких двух разных решений быть не может. В терминах теоремы Кронекера — Капелли это тот самый случай, когда ранг матрицы равен рангу расширенной матрицы и равен числу неизвестных: система совместна и определённа.

Почему именно слева, и что сломается справа

Это место стоит разобрать медленно, потому что оно даёт очень частую ошибку.

Попробуем домножить $Ax = b$ на $A^{-1}$ справа. Получится выражение

и здесь начинается арифметика размеров. Матрица $A$ имеет размер $n \times n$, вектор $x$ — размер $n \times 1$, значит их произведение $Ax$ имеет размер $n \times 1$. Теперь мы хотим умножить объект размера $n \times 1$ на матрицу размера $n \times n$. Правило умножения матриц требует, чтобы число столбцов левого множителя совпадало с числом строк правого: слева столбцов $1$, справа строк $n$. При $n > 1$ это не определено вообще — такого произведения не существует. То же самое с правой частью: $bA^{-1}$ — это $(n\times1)\cdot(n\times n)$, тоже несуществующая операция.

То есть домножение справа на $A^{-1}$ в системе $Ax=b$ ломается не «даёт неправильный ответ», а «не является математической записью». Это хорошая новость: ошибку сразу видно по размерам, если привыкнуть их проверять.

Но есть случай, где размеры сходятся и ошибка становится незаметной, — матричное уравнение, где неизвестное само квадратная матрица. Если $AX = B$, где все три матрицы имеют размер $n\times n$, то запись $XA^{-1}$ и $BA^{-1}$ вполне определена. И вот здесь домножение не с той стороны даёт формально корректное, но неверное выражение: из $AX = B$ домножением справа получается $AXA^{-1} = BA^{-1}$, а $AXA^{-1}$ никак не упрощается до $X$ — сократить $A$ слева и $A^{-1}$ справа невозможно, потому что между ними стоит $X$, а переставлять множители нельзя.

Простое правило, которое стоит запомнить: множитель убирают с той стороны, с которой он стоит. В $AX = B$ матрица $A$ стоит слева от неизвестного — домножаем слева, получаем $X = A^{-1}B$. В $XA = B$ матрица $A$ стоит справа — домножаем справа, получаем $X = BA^{-1}$. Система $Ax = b$ — это частный случай первого варианта, где $X$ состоит из одного столбца.

Разбор примеров

Пример 1: система 3×3 с целочисленной обратной

Решим ту самую систему, с которой начали:

Шаг 1. Проверяем применимость метода. Считаем определитель разложением по первой строке:

$\det A = -1 \ne 0$ — матрица невырождена, метод применим, решение существует и единственно.

Шаг 2. Находим $A^{-1}$. Пользуемся уже знакомой техникой из урока про обратную матрицу; здесь нас интересует только результат:

Проверим, что это действительно обратная, — умножим на $A$ и убедимся, что получится $E$. Первая строка $A^{-1}$ на первый столбец $A$: $4\cdot2 + 3\cdot(-3) + (-1)\cdot(-2) = 8 - 9 + 2 = 1$ ✅. Первая строка $A^{-1}$ на второй столбец $A$: $4\cdot1 + 3\cdot(-1) + (-1)\cdot 1 = 4 - 3 - 1 = 0$ ✅. Остальные семь произведений считаются так же и дают единичную матрицу.

Шаг 3. Умножаем $A^{-1}$ на $b$. Это обычное умножение матрицы на вектор — каждая координата ответа есть скалярное произведение строки $A^{-1}$ на $b$:

Шаг 4. Проверка подстановкой в исходную систему — обязательный шаг, который занимает двадцать секунд и ловит почти любую арифметическую ошибку:

Ответ: $x_1 = 2$, $x_2 = 3$, $x_3 = -1$.

Пример 2: система 3×3 с дробной обратной

Шаг 1. Определитель.

Метод применим.

Шаг 2. Обратная матрица.

Шаг 3. Умножаем на $b = (5,\ 0,\ 6)^T$. Удобно не расписывать дроби, а вынести $\frac15$ за скобку и работать с целыми числами:

Шаг 4. Проверка:

Ответ: $x = (2,\ 1,\ -1)^T$.

Обрати внимание на технический приём из шага 3: вместо того чтобы делить каждый элемент союзной матрицы на определитель и потом возиться с дробями, множитель $\frac{1}{\det A}$ выносится наружу, всё умножение делается в целых числах, а деление выполняется один раз в самом конце. На матрицах $3\times3$ и больше это экономит массу арифметических ошибок.

Пример 3: когда метод не работает

Шаг 1. Определитель. Вторая строка матрицы коэффициентов — это первая, умноженная на 2. По свойству определителей, если одна строка пропорциональна другой, определитель равен нулю:

Шаг 2. Останавливаемся. Матрица вырождена, $A^{-1}$ не существует, и формула $x = A^{-1}b$ не имеет смысла. Дальше считать нечего — никакие ухищрения не позволят применить матричный метод к этой системе.

Это не значит, что у системы нет решений. Здесь второе уравнение — точная копия первого, умноженная на 2 (включая правую часть: $8 = 2\cdot 4$), поэтому оно не несёт новой информации, и система сводится к двум независимым уравнениям с тремя неизвестными — у неё бесконечно много решений. Но матричный метод их не найдёт: он в принципе умеет находить только единственное решение и только у квадратных невырожденных систем. Для вырожденного случая нужен метод Гаусса, а если правая часть нулевая — техника из следующего урока про однородные системы.

Вывод из примера: проверка $\det A \ne 0$ — это не формальность в начале решения, а собственно критерий применимости метода. Если определитель равен нулю, метод не «даёт плохой ответ», он не применим вовсе.

Почему это важно

Матричный метод — это в первую очередь способ думать, а не способ считать. Формула $x = A^{-1}b$ говорит: решение системы линейно зависит от правой части. Удвоил $b$ — удвоилось решение. Сложил две правые части — сложились решения. Эта линейность неочевидна, если смотреть на систему как на набор уравнений, но становится тривиальной, как только ты видишь оператор $A^{-1}$, действующий на $b$.

Из той же формулы сразу читается и главный вопрос устойчивости, которому посвящена вторая половина урока. Если $b$ известна неточно — а в реальных задачах она всегда известна неточно, это же измерения, — то ошибка $\delta b$ превращается в ошибку решения $\delta x = A^{-1}\,\delta b$. То есть матрица $A^{-1}$ работает усилителем ошибки, и вопрос «насколько можно верить ответу» превращается в вопрос «насколько сильно $A^{-1}$ растягивает векторы». Ровно это мы и научимся измерять.

Системы с несколькими правыми частями: $AX = B$

Интуиция: одна фабрика, много сценариев

Представь цех, который выпускает три вида изделия. На каждое изделие уходит фиксированное количество трёх ресурсов — металла, пластика и человеко-часов. Расход не меняется: технология одна. А вот запасы ресурсов на складе меняются каждую неделю, и каждую неделю нужно заново отвечать на вопрос «какой план выпуска полностью выберет склад?».

Математически это система $Ax = b$, где матрица $A$ (расход ресурсов на единицу изделия) всегда одна и та же, а вектор $b$ (запасы) каждый раз новый. За год таких вопросов набирается пятьдесят два.

Решать пятьдесят два раза с нуля методом Гаусса — значит пятьдесят два раза проделать одну и ту же работу по приведению одной и той же матрицы к ступенчатому виду. Вся тяжёлая часть вычислений зависит только от $A$ и вообще не зависит от $b$ — и её обидно повторять.

Вот здесь у матричного метода появляется настоящее преимущество: обратную матрицу считают один раз, а дальше каждый новый ответ — это одно умножение матрицы на вектор, то есть буквально $n^2$ умножений и сложений. Для $n = 3$ это девять умножений вместо полного прогона Гаусса.

Постановка

Пусть у нас есть $k$ систем с одной и той же матрицей $A$ и разными правыми частями $b^{(1)}, b^{(2)}, \dots, b^{(k)}$. Соберём правые части в матрицу $B$ размера $n \times k$, поставив их столбцами, а искомые решения — в матрицу $X$ того же размера:

Тогда все $k$ систем записываются одним равенством:

Это верно потому, что умножение матриц действует на столбцы независимо: $j$-й столбец произведения $AX$ равен $A$, умноженной на $j$-й столбец $X$. То есть равенство $AX = B$ по столбцам распадается ровно на $k$ отдельных систем $A x^{(j)} = b^{(j)}$.

Решение получается домножением слева (сторона та же, что и раньше, и по той же причине):

Обрати внимание на ракурс. В уроке про обратную матрицу уравнение $AX = B$ было самостоятельным объектом — «найди неизвестную матрицу». Здесь это то же уравнение, но читается оно иначе: не «одно матричное уравнение», а «пачка обычных систем, упакованная в одну запись». Столбцы $B$ — это независимые сценарии, столбцы $X$ — независимые ответы. Ни один столбец не влияет на другой.

Разбор примера: три сценария поставок

Пусть цех выпускает три изделия. Расход трёх ресурсов на единицу каждого изделия задан матрицей

Строка — ресурс, столбец — изделие. Скажем, элемент $a_{23} = 3$ означает: на одну единицу третьего изделия уходит 3 единицы второго ресурса.

На три ближайшие недели склад обещает такие запасы (столбцы — недели):

Нужно найти планы выпуска на каждую из трёх недель.

Шаг 1. Определитель. Разложим по первой строке:

$\det A = 1 \ne 0$ — метод применим. (Определитель, равный единице, — приятная случайность выбранного примера: обратная матрица получится целочисленной.)

Шаг 2. Обратная матрица — считаем её один раз для всех трёх недель:

Быстрая проверка: первая строка $A^{-1}$ на первый столбец $A$ даёт $3\cdot1 + (-3)\cdot1 + 1\cdot1 = 1$ ✅, первая строка на второй столбец: $3\cdot1 + (-3)\cdot2 + 1\cdot3 = 3 - 6 + 3 = 0$ ✅.

Шаг 3. Одно матричное умножение вместо трёх решений систем:

Считаем первый столбец результата:

Второй столбец:

Третий столбец:

Шаг 4. Проверка. Подставим первый столбец в исходную систему: $20 + 30 + 50 = 100$ ✅, $20 + 2\cdot30 + 3\cdot50 = 20 + 60 + 150 = 230$ ✅, $20 + 3\cdot30 + 6\cdot50 = 20 + 90 + 300 = 410$ ✅. Второй и третий столбцы проверяются так же.

Ответ: на первой неделе выпускаем $(20,\ 30,\ 50)$, на второй $(50,\ 10,\ 20)$, на третьей $(10,\ 40,\ 30)$ единиц трёх изделий соответственно.

Когда $A^{-1}$ действительно окупается — и когда нет

Давай честно посчитаем. Мы будем измерять работу в арифметических операциях (умножениях и сложениях); для больших $n$ это адекватная модель времени работы.

• Метод Гаусса для одной системы: прямой ход стоит примерно $\dfrac{2n^3}{3}$ операций, обратный ход — ещё около $2n^2$. Итого $\approx \dfrac{2n^3}{3}$.

• Решать $k$ систем методом Гаусса «с нуля» каждый раз: $k \cdot \dfrac{2n^3}{3}$. Кубическая работа повторяется $k$ раз.

• Построить $A^{-1}$ явно (методом Гаусса — Жордана над расширенной матрицей $(A \mid E)$) стоит примерно $2n^3$ операций. Дальше каждое решение — это $A^{-1}b^{(j)}$, то есть $2n^2$ операций. Итого $2n^3 + k\cdot 2n^2$.

Подставим $n = 1000$ и посмотрим на числа (операций, порядок величины):

Картина ясная. При одной правой части явная обратная матрица проигрывает Гауссу примерно втрое: ты делаешь втрое больше работы, чтобы получить тот же ответ. Но начиная примерно с $k \approx 3$ ситуация переворачивается, а при $k = 100$ явная $A^{-1}$ выигрывает уже в тридцать раз. Это и есть честный аргумент «за» матричный метод: он окупается, когда матрица одна, а правых частей много.

Теперь честный аргумент «против», без которого картина была бы неполной. Идея «тяжёлую работу над $A$ сделать один раз» правильная, но явная обратная матрица — не лучший способ её реализовать. Существует LU-разложение: матрица $A$ один раз раскладывается в произведение двух треугольных множителей за те же $\frac{2n^3}{3}$ операций, а дальше каждая новая правая часть обрабатывается двумя подстановками за $2n^2$ операций. Итого $\frac{2n^3}{3} + k\cdot 2n^2$ — это дешевле явной обратной матрицы при любом $k$, и при этом численно точнее. Технику LU мы разбирать не будем (это тема более поздних курсов), но знать про её существование нужно: именно она стоит внутри numpy.linalg.solve, и именно поэтому в NumPy можно передать в solve сразу матрицу правых частей:

import numpy as np

A = np.array([[1., 1., 1.],
              [1., 2., 3.],
              [1., 3., 6.]])
B = np.array([[100., 80., 80.],
              [230., 130., 180.],
              [410., 200., 310.]])

X = np.linalg.solve(A, B)   # сразу три системы, одно разложение
print(X)
# [[20. 50. 10.]
#  [30. 10. 40.]
#  [50. 20. 30.]]

Один вызов solve с матрицей $B$ раскладывает $A$ один раз и решает все три системы — то есть даёт ровно ту экономию, ради которой мы затевали явную $A^{-1}$, но дешевле и устойчивее.

Остаётся один сценарий, где явная обратная матрица нужна именно как объект, а не как способ решить систему: когда сами её элементы что-то значат. Классический пример — модель «затраты — выпуск» Василия Леонтьева в экономике, где матрица $(E - A)^{-1}$ называется матрицей полных затрат, и её элемент $c_{ij}$ прямо читается как «сколько продукции отрасли $i$ нужно суммарно, прямо и косвенно, чтобы выдать единицу конечного продукта отрасли $j$». Там обратную матрицу считают не чтобы решить систему, а чтобы посмотреть на её числа. Похожая история в статистике: диагональ матрицы $(X^TX)^{-1}$ даёт дисперсии оценок коэффициентов регрессии, и её иногда действительно нужно вычислить целиком.

Почему это важно

В машинном обучении сценарий «одна матрица, много правых частей» встречается постоянно, причём часто в неявном виде. Обучение модели гребневой регрессии на сетке из двадцати значений параметра регуляризации — это двадцать систем, но, к сожалению, с разными матрицами. А вот вычисление предсказаний гауссовского процесса в тысяче новых точек — это уже честная тысяча правых частей при одной и той же ковариационной матрице. Решение нескольких задач оптимального управления с одной моделью динамики — то же самое. В таких задачах разница между «переразложить матрицу тысячу раз» и «разложить один раз» — это разница между часом и секундой.

И ещё один момент, который стоит держать в голове: $AX = B$ с $B = E$ — это в точности задача нахождения обратной матрицы, $X = A^{-1}E = A^{-1}$. То есть обращение матрицы — частный случай решения системы с $n$ правыми частями, а не наоборот. Это переворачивает привычную иерархию: не «умею обращать, значит умею решать», а «умею решать, значит умею обращать, если очень надо».

Три метода в очной ставке: Гаусс, Крамер, матричный

К этому моменту у тебя в руках три способа решить систему линейных уравнений. Пора разложить их по полкам честно, с числами, а не лозунгами.

Область применимости

Это первое, что разводит методы, и разводит жёстко.

Метод Гаусса работает с любой системой: квадратной и прямоугольной, совместной и несовместной, определённой и неопределённой. Он не требует, чтобы решение существовало, — он просто приводит систему к ступенчатому виду и честно сообщает, что получилось: единственное решение, бесконечное семейство с параметрами или противоречие вида $0 = 1$.

Правило Крамера требует, чтобы матрица была квадратной и невырожденной ($\Delta \ne 0$). При $\Delta = 0$ оно не даёт ответа, а только сигнализирует, что случай особый.

Матричный метод требует ровно того же: квадратная матрица, $\det A \ne 0$. Область применимости у Крамера и матричного метода совпадает буквально — это два разных способа записать одно и то же решение. Более того, они математически тождественны: если расписать $x = \frac{1}{\det A}\operatorname{adj}(A)\, b$ покоординатно, каждая координата окажется в точности крамеровской дробью $\Delta_i/\Delta$. Разные внешне, один и тот же объект внутри.

Стоимость

Считаем арифметические операции для системы $n\times n$. Формулы такие:

• Гаусс: $\dfrac{2n^3}{3} + 2n^2$;

• Матричный метод (явная $A^{-1}$ через Гаусса — Жордана, затем $A^{-1}b$): $2n^3 + 2n^2$;

• Крамер, если каждый определитель считать методом Гаусса: $(n+1)\cdot\dfrac{2n^3}{3}$;

• Крамер, если каждый определитель раскладывать по строке (как учат руками): порядка $(n+1)\cdot n!$.

Числа:

Читаем таблицу. На $n = 3$ разница между методами несущественна — три десятка операций против семи десятков, для ручного счёта роли не играет, выбирай что удобнее. На $n = 10$ матричный метод уже в 2,5 раза дороже Гаусса, а Крамер «как учат руками» — в сорок пять тысяч раз дороже. На $n = 20$ наивный Крамер требует $5\cdot10^{19}$ операций: на машине, делающей миллиард операций в секунду, это 1622 года. Гауссу на ту же систему нужно 6133 операции — примерно шесть микросекунд.

Отдельно про матричный метод и множество правых частей: если правых частей $k$, его стоимость $2n^3 + k\cdot 2n^2$ против $k\cdot\frac{2n^3}{3}$ у наивного повторного Гаусса. Точка окупаемости — около $k = 3$, дальше преимущество растёт (числа мы посчитали в предыдущем разделе).

Устойчивость

Это измерение, которое в учебниках упоминают реже всего, а на практике оно решает.

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

Матричный метод уступает по двум причинам сразу. Во-первых, он делает больше арифметических операций, а каждая операция в машинной арифметике вносит крошечную ошибку округления, и их больше накапливается. Во-вторых, и это важнее, он проходит через промежуточный объект $A^{-1}$, элементы которого сами вычислены с ошибкой; затем эта уже испорченная матрица умножается на $b$, и ошибка умножается вместе с ней. Метод Гаусса такого промежуточного объекта не создаёт вовсе. Насколько это существенно — мы сейчас измерим экспериментально.

Правило Крамера численно хуже всех. Каждая координата ответа получается делением одного определителя на другой, а определители — величины, крайне чувствительные к округлению: они масштабируются как $n$-я степень элементов и легко переполняются или обнуляются. Для $n$ больше 3–4 правило Крамера как вычислительный метод не используют нигде.

Сводная таблица

Почему это важно

Главный вывод таблицы не «Гаусс лучше всех», а более тонкий: у методов разные задачи. Правило Крамера — инструмент вывода формул: когда нужно понять, как решение зависит от параметра, крамеровская дробь показывает это явно, а прогон Гаусса — нет. Матричный метод — инструмент рассуждения: он делает видимой линейность решения по правой части и превращает систему в «применение оператора $A^{-1}$». Метод Гаусса — инструмент счёта: он универсален, дёшев и устойчив, и именно он живёт внутри библиотек.

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

Насколько можно доверять ответу: нормы и число обусловленности

Интуиция: матрица как усилитель ошибки

Все примеры выше были стерильными: правая часть задана точными целыми числами, определитель считается без округлений, ответ выходит красивым. В реальной задаче не так. Правая часть $b$ — это измерения: показания датчиков, суммы по чекам, замеры расхода. Они известны с погрешностью. Пусть погрешность крошечная, скажем 0,01 %. Вопрос: какой будет погрешность в ответе $x$?

Наивный ответ — «тоже около 0,01 %». Он неверен, и вот почему. Из формулы $x = A^{-1}b$ мгновенно следует, что если правая часть изменилась на $\delta b$, то решение изменилось на

Матрица $A^{-1}$ действует на вектор ошибки так же, как на любой другой вектор: растягивает его. Если $A^{-1}$ умеет растягивать векторы в десять тысяч раз, то ошибка в 0,01 % на входе превратится в ошибку в 100 % на выходе. И это не патология, а обычное поведение: сейчас мы увидим матрицу $2\times2$ из почти целых чисел, которая делает ровно это.

Чтобы говорить о «растяжении» количественно, нужна мера длины вектора и мера «силы растяжения» матрицы. Это и есть нормы.

Норма вектора: сколько это «длинно»

Определение: Норма вектора — это число $\|x\| \ge 0$, играющее роль его длины. Основные три:

$$\|x\|_2 = \sqrt{x_1^2 + x_2^2 + \dots + x_n^2} \qquad \text{(евклидова, «обычная» длина)}$$$$\|x\|_1 = |x_1| + |x_2| + \dots + |x_n| \qquad \text{(«манхэттенская»)}$$$$\|x\|_\infty = \max_i |x_i| \qquad \text{(максимальная координата по модулю)}$$

Основная для нас — евклидова $\|x\|_2$: это ровно та длина вектора, которую ты считал в уроках про векторы, просто записанная в $n$ измерениях. Например, для $x = (3,\ -4,\ 12)$:

Все три нормы дают разные числа, но ведут себя одинаково в главном: они равны нулю только у нулевого вектора, умножаются на $|\alpha|$ при умножении вектора на $\alpha$ и удовлетворяют неравенству треугольника $\|u + v\| \le \|u\| + \|v\|$. Для наших целей выбор нормы почти не важен — важно, что мы можем сравнить «размер» ошибки с «размером» самого вектора.

Ключевая величина — не абсолютная ошибка, а относительная:

Именно относительная погрешность отвечает на человеческий вопрос «сколько верных знаков». Относительная погрешность $10^{-3}$ означает примерно три верных десятичных знака, $10^{-6}$ — шесть, и так далее.

Норма матрицы: максимальное растяжение

Теперь нужна мера того, насколько сильно матрица растягивает векторы. Матрица $A$ разные векторы растягивает по-разному: какой-то удлиняет вдвое, какой-то, наоборот, укорачивает. Разумно взять худший случай.

Определение: Нормой матрицы $A$ (согласованной с векторной нормой $\|\cdot\|$) называется максимальный коэффициент растяжения:

$$\|A\| = \max_{x \ne 0} \frac{\|Ax\|}{\|x\|}$$

То есть $\|A\|$ — это ответ на вопрос: «во сколько раз максимум может удлиниться вектор, если применить к нему $A$?»

Из определения сразу следует главное рабочее неравенство, которым мы будем пользоваться постоянно:

и его аналог для произведения матриц: $\|AB\| \le \|A\|\cdot\|B\|$.

Считать максимум по всем векторам вручную неудобно, поэтому есть готовые формулы. Две из них считаются устно:

Норма $\|A\|_\infty$ согласована с векторной нормой $\|x\|_\infty$, норма $\|A\|_1$ — с $\|x\|_1$. А вот $\|A\|_2$ (согласованная с обычной евклидовой длиной) вручную не считается: она равна наибольшему сингулярному числу матрицы, и это тема будущих уроков. Практический вывод: руками считаем через $\|\cdot\|_\infty$, а если нужно точное значение в $\|\cdot\|_2$, спрашиваем у numpy.linalg.norm(A, 2).

Пример. Для $A = \begin{pmatrix} 2 & -3 \\ 4 & 1 \end{pmatrix}$:

суммы модулей по строкам: $|2| + |-3| = 5$ и $|4| + |1| = 5$, значит $\|A\|_\infty = 5$;

суммы модулей по столбцам: $|2| + |4| = 6$ и $|-3| + |1| = 4$, значит $\|A\|_1 = 6$.

(Для сравнения, numpy даёт $\|A\|_2 \approx 4{,}515$ — как и полагается, значение между ними по порядку, но не совпадающее ни с одним.)

Число обусловленности

Теперь собираем всё вместе. Пусть $Ax = b$ и правая часть возмущена: вместо $b$ мы решаем систему с $b + \delta b$ и получаем $x + \delta x$. Вычитая одно равенство из другого, получаем $A\,\delta x = \delta b$, то есть $\delta x = A^{-1}\delta b$. Оценим:

С другой стороны, из $b = Ax$ следует $\|b\| \le \|A\|\cdot\|x\|$, то есть $\dfrac{1}{\|x\|} \le \dfrac{\|A\|}{\|b\|}$. Перемножаем два неравенства:

Множитель, который выскочил перед относительной погрешностью входа, и есть главный герой урока.

Определение: Числом обусловленности невырожденной матрицы $A$ называется

$$\operatorname{cond}(A) = \|A\|\cdot\|A^{-1}\|$$

Оно даёт оценку усиления относительной погрешности:

$$\frac{\|\delta x\|}{\|x\|} \ \le\ \operatorname{cond}(A)\cdot\frac{\|\delta b\|}{\|b\|}$$

Читается это так: «ошибка в правой части может вырасти в $\operatorname{cond}(A)$ раз, но не больше». Матрицы с маленьким $\operatorname{cond}$ называют хорошо обусловленными, с большим — плохо обусловленными.

Три свойства, которые полезно понимать сразу.

Первое: $\operatorname{cond}(A) \ge 1$ всегда. Действительно, $1 = \|E\| = \|A A^{-1}\| \le \|A\|\cdot\|A^{-1}\| = \operatorname{cond}(A)$. Идеальный минимум $\operatorname{cond} = 1$ достигается, например, у единичной матрицы и у матриц поворота: они вообще не искажают длины, и ошибка на выходе ровно такая же, как на входе.

Второе: $\operatorname{cond}$ не зависит от масштаба. Если умножить всю матрицу на число $\alpha \ne 0$, то $\|\alpha A\| = |\alpha|\,\|A\|$, а $\|(\alpha A)^{-1}\| = \frac{1}{|\alpha|}\|A^{-1}\|$, и множители сокращаются: $\operatorname{cond}(\alpha A) = \operatorname{cond}(A)$. Это правильно: умножение всех уравнений на 1000 не делает систему ни проще, ни сложнее.

Третье: большое $\operatorname{cond}$ означает близость к вырожденности. Если $\det A$ близок к нулю, то элементы $A^{-1}$ (в которых стоит деление на $\det A$) огромны, значит $\|A^{-1}\|$ огромна, значит $\operatorname{cond}$ огромно. При этом связь между $\det A$ и обусловленностью не прямая: сам по себе маленький определитель ничего не доказывает (у матрицы $0{,}001\cdot E$ определитель $10^{-9}$ при $n=3$, а $\operatorname{cond} = 1$ — она идеальна). Именно поэтому определитель — плохой индикатор численных проблем, а $\operatorname{cond}$ — хороший.

Есть и удобное «правило большого пальца» для оценки потери точности:

Числа с плавающей точкой двойной точности (тип float64) хранят около 16 верных десятичных знаков. Значит, при $\operatorname{cond} = 10^{4}$ в ответе останется примерно $16 - 4 = 12$ верных знаков — не страшно. При $\operatorname{cond} = 10^{12}$ останется четыре. При $\operatorname{cond} = 10^{16}$ не останется ни одного, и любой ответ, который выдаст компьютер, будет случайным шумом.

Забегая вперёд одной фразой: у нормы $\|\cdot\|_2$ число обусловленности равно отношению наибольшего сингулярного числа матрицы к наименьшему, $\operatorname{cond}_2(A) = \sigma_{\max}/\sigma_{\min}$, — это самое употребительное определение в численных библиотеках, но сингулярные числа появятся у тебя в курсе позже, а через нормы всё считается уже сейчас.

Разбор примеров

Пример 1: система $2\times2$, где ответ улетает

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

Решение. Вычтем первое уравнение из второго: $0{,}0001\,x_2 = 0{,}0001$, откуда $x_2 = 1$ и $x_1 = 1$. Решение $x = (1,\ 1)$, ровное и симпатичное.

Теперь испортим правую часть. Пусть в результате измерения второй свободный член оказался не $2{,}0001$, а $2{,}0002$ — разница в одну десятитысячную:

Вычитаем так же: $0{,}0001\,x_2 = 0{,}0002$, откуда $x_2 = 2$ и $x_1 = 0$.

Останови взгляд на этом. Решение было $(1,\ 1)$, стало $(0,\ 2)$. Ни одна координата не сохранилась даже приблизительно.

Посчитаем усиление в числах. В евклидовой норме:

Фактическое усиление: $\dfrac{1}{3{,}54\cdot10^{-5}} \approx 28\,285$ раз.

Проверим, предсказывает ли это число обусловленности. Матрица и её обратная:

Считаем нормы по строкам:

numpy.linalg.cond(A) в евклидовой норме даёт $40\,002{,}0$ — практически то же самое.

Фактическое усиление $28\,285$ оказалось меньше оценки $40\,002$ — так и должно быть, ведь $\operatorname{cond}$ даёт верхнюю границу для худшего направления ошибки. Проверим, достижима ли граница: если сдвинуть $b$ ровно на 0,01 % от её длины, но в самом «неудачном» направлении, численный эксперимент даёт новое решение $x' = (-3{,}0003;\ 5{,}0001)$, то есть относительное изменение решения $4{,}0002$, или 400,02 %. Усиление: $4{,}0002 / 10^{-4} = 40\,002$ — граница достигнута с точностью до последнего знака.

Вот и весь смысл понятия: сотая доля процента на входе, четыреста процентов на выходе.

Откуда взялась беда? Геометрически два уравнения задают две прямые на плоскости, а решение — точку их пересечения. Здесь прямые почти параллельны: их коэффициенты $(1;1)$ и $(1;1{,}0001)$ различаются в четвёртом знаке. Точка пересечения почти параллельных прямых чудовищно чувствительна к их положению: сдвинь одну прямую на волосок — и точка пересечения уедет на километр. Плохая обусловленность — это алгебраическое имя для «почти параллельных уравнений».

Пример 2: матрица Гильберта $4\times4$

Классический тестовый случай. Матрица Гильберта — это $h_{ij} = \dfrac{1}{i + j - 1}$:

Ничего подозрительного: симметричная, все элементы положительные, порядка единицы. Определитель $\det H = \dfrac{1}{6\,048\,000} \approx 1{,}65\cdot10^{-7}$ — маленький, но не нулевой. Обратная матрица целочисленная и огромная:

(В евклидовой норме numpy даёт $\operatorname{cond}_2(H) = 15\,514$ — тот же порядок.)

Ставим эксперимент. Пусть точное решение системы $Hx = b$ — это $x = (1,\ 1,\ 1,\ 1)$. Тогда правая часть равна суммам строк:

Теперь представим самую обычную ситуацию: правую часть кто-то записал с тремя знаками после запятой. Это ведь очень мелкое огрубление:

Относительное изменение правой части: $\dfrac{\|\tilde b - b\|_2}{\|b\|_2} = 2{,}45\cdot10^{-4}$, то есть 0,0245 %.

Решаем систему с этой округлённой правой частью (расчёт в numpy):

Вместо $(1;\ 1;\ 1;\ 1)$. Третья координата даже сменила знак. Относительное изменение решения — $75{,}2\ \%$, фактическое усиление $3065$ раз (в пределах оценки $\operatorname{cond}_2 = 15\,514$, как и должно быть).

А если округлить правую часть до двух знаков — $\tilde b = (2{,}08;\ 1{,}28;\ 0{,}95;\ 0{,}76)$, изменение на 0,17 %, — решение станет

От исходного $(1;1;1;1)$ не осталось ничего.

И честная оговорка, которая важнее самого эффекта. Не всякое возмущение раздувается. Если умножить всю правую часть на $1{,}0001$ (то есть изменить $b$ на 0,01 %, но строго вдоль самого вектора $b$), решение станет ровно $(1{,}0001;\ 1{,}0001;\ 1{,}0001;\ 1{,}0001)$ — усиление в точности единица, никакой катастрофы. Число обусловленности описывает худший случай, а не типичный. Проблема в том, что реальная ошибка измерений — это, как правило, произвольное направление, и рано или поздно она попадёт в опасное.

Пример 3: хорошо обусловленная система для контраста

Чтобы не создалось впечатления, будто всё плохо всегда, посмотрим на нормальную матрицу:

Она симметрична, с явным преобладанием диагонали. $\det A = 97$, обратная:

Три с небольшим. Это значит, что относительная погрешность решения не превысит утроенной погрешности входных данных: измерил правую часть с точностью 1 % — получишь ответ с точностью не хуже 3,2 %. С такой матрицей можно работать спокойно, и никакие численные тонкости здесь роли не играют.

Диагональное преобладание (когда $|a_{ii}|$ больше суммы модулей остальных элементов строки) — хороший практический признак: такие матрицы почти всегда хорошо обусловлены. В плохо обусловленных матрицах, наоборот, строки почти пропорциональны друг другу.

Почему это важно

Число обусловленности разделяет два принципиально разных вида проблем, которые новички постоянно путают.

Первый вид — проблема алгоритма: ты выбрал плохой метод, и он накопил лишнюю ошибку. Это лечится сменой метода (например, переходом с явной обратной матрицы на solve).

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

Различать эти два случая — практический навык. Если у тебя $\operatorname{cond}(A) = 10^{14}$, менять inv на solve бесполезно: надо менять постановку задачи — регуляризовать, отбросить лишние признаки, перемасштабировать переменные. Если $\operatorname{cond}(A) = 20$, а ответы всё равно скачут, то виноват не «плохой» матан, а твой код.

Почему в коде никогда не пишут inv(A) @ b

Интуиция: лишний посредник

Две строчки кода:

x = np.linalg.solve(A, b)      # так правильно
x = np.linalg.inv(A) @ b       # так не надо

Математически они означают одно и то же. Вычислительно — нет, и разница не косметическая.

Аналогия. Тебе нужно узнать, сколько будет $\dfrac{347}{53}$. Первый способ: поделить и получить $6{,}5472$. Второй способ: сначала вычислить $\dfrac{1}{53} = 0{,}018868$, округлив до шести знаков, а потом умножить $347 \cdot 0{,}018868 = 6{,}54720$. Ответы почти совпадут, но во втором способе ты сделал лишнюю операцию округления и потратил больше действий — ради ровно того же результата. С матрицами всё то же самое, только «лишнее округление» происходит $n^2$ раз (по разу на каждый элемент $A^{-1}$), а «лишняя работа» составляет две трети всех вычислений.

Разберём три причины по отдельности и измерим каждую.

Причина 1: точность

Вот эксперимент. Берём матрицы Гильберта разного размера, задаём точное решение $x = (1, 1, \dots, 1)$, вычисляем $b = Hx$, а затем восстанавливаем $x$ двумя способами и смотрим на относительную ошибку $\dfrac{\|\hat x - x\|_2}{\|x\|_2}$:

Читаем. Во-первых, видно, как работает правило потери знаков: при $n = 10$ и $\operatorname{cond} \approx 10^{13}$ мы теряем 13 знаков из 16, остаётся три — и ошибка solve действительно порядка $10^{-4}$. Всё сходится. Во-вторых, inv стабильно даёт ошибку в 15–70 раз больше при абсолютно одинаковых входных данных. Это чистая потеря точности на пустом месте: те же данные, тот же ответ, просто хуже посчитано.

Ещё нагляднее ситуация с невязкой — величиной $\dfrac{\|A\hat x - b\|}{\|b\|}$, которая показывает, насколько найденный ответ вообще удовлетворяет исходным уравнениям. Возьмём случайную матрицу $200\times200$ с $\operatorname{cond} = 10^{10}$:

• ошибка решения: solve даёт $9{,}0\cdot10^{-8}$, inv @ b даёт $8{,}5\cdot10^{-7}$ (в 9,5 раза хуже);

• невязка: solve даёт $1{,}0\cdot10^{-15}$, inv @ b даёт $1{,}3\cdot10^{-7}$ — в сто миллионов раз хуже.

Разрыв в невязке огромен, и он показателен. solve устроен так, что найденный им $\hat x$ практически идеально удовлетворяет системе — невязка на уровне машинного эпсилона. Это гарантия обратной устойчивости: даже если ответ далёк от истинного (потому что задача плохо обусловлена), он является точным решением почти той же системы. У inv @ b этой гарантии нет: полученный вектор не удовлетворяет исходным уравнениям даже приблизительно. То есть solve терпит поражение только там, где задача сама безнадёжна, а inv подводит дополнительно, от себя.

Причина 2: скорость

Теоретически мы уже посчитали: явная обратная матрица стоит $2n^3$ операций против $\frac{2n^3}{3}$ у Гаусса, то есть должна быть примерно втрое медленнее. Проверим на практике. Замер на случайных матрицах с диагональным преобладанием, numpy поверх стандартного LAPACK, время одного вызова (минимум из серии повторов):

Теория предсказывала трёхкратное отставание — практика дала от 2,6 до 3,8 раз. Совпадение хорошее (разброс объясняется тем, что современные библиотеки по-разному распараллеливают операции разного типа).

Три раза — это не «мелочь, о которой не стоит думать». Если решение системы стоит внутри цикла обучения и вызывается миллион раз, разница между 11 мс и 32 мс — это разница между тремя часами и девятью.

Причина 3: память

Пункт, о котором вспоминают реже всего, а он бывает решающим.

Матрица $A$ размера $n \times n$ в реальных задачах часто разреженная: почти все элементы — нули. Матрица конечных элементов, матрица смежности графа, матрица ковариаций с обрезанными связями — во всех этих случаях ненулевых элементов порядка $n$, а не $n^2$. Разреженную матрицу $10^6 \times 10^6$ с десятью ненулями в строке можно спокойно держать в памяти: это десять миллионов чисел, около 80 мегабайт.

А вот обратная к разреженной матрице почти всегда плотная — нули в ней не сохраняются. Обратная к той же матрице потребовала бы $10^{12}$ чисел, то есть восемь терабайт. Она физически не помещается никуда. При этом решить систему с разреженной матрицей вполне реально — специализированные разреженные решатели (scipy.sparse.linalg.spsolve) это делают, никогда не строя обратную.

То же самое, но мягче, действует и на плотных матрицах: solve работает «на месте», а inv требует дополнительный массив на $n^2$ чисел.

Итоговое правило

Правило solve, а не inv: если тебе нужен вектор $x$, решающий $Ax = b$, вызывай функцию решения системы, а не строй обратную матрицу. Обратная матрица нужна только тогда, когда тебе нужны её собственные элементы как самостоятельные числа.

x = np.linalg.solve(A, b)          # ✅
x = np.linalg.inv(A) @ b           # ❌
X = np.linalg.solve(A, B)          # ✅ несколько правых частей сразу
X = np.linalg.inv(A) @ B           # ❌

Это правило прямо записано в документации NumPy и MATLAB, в книге Голуба и Ван Лоуна «Matrix Computations» и в кодстайлах большинства ML-команд. Теперь ты знаешь не только формулировку, но и три причины за ней.

Вырожденный и почти вырожденный случай

Отдельный сюжет — что происходит, когда матрица не обратима вовсе или почти не обратима. Разница между этими двумя случаями в коде принципиальная, и она контринтуитивна.

Честно вырожденная матрица падает с ошибкой, и это хорошо. Возьмём

Вторая строка ровно вдвое больше первой, $\det A = 0$. Что скажет код:

A = np.array([[1., 2.], [2., 4.]])
b = np.array([3., 6.])
np.linalg.solve(A, b)
# numpy.linalg.LinAlgError: Singular matrix

Исключение LinAlgError: Singular matrix — это лучший исход из возможных. Программа остановилась, ты немедленно узнал о проблеме, посмотрел на матрицу и обнаружил, что второе уравнение дублирует первое. У этой системы, кстати, бесконечно много решений (все точки прямой $x_1 + 2x_2 = 3$), но solve не умеет их выдавать — он для однозначно разрешимых систем.

Почти вырожденная матрица не падает — и это плохо. Испортим предыдущую матрицу на одну стомиллионную:

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

A = np.array([[1., 2.], [2., 4.0000001]])
np.linalg.solve(A, np.array([3., 6.]))          # -> [3., 0.]
np.linalg.solve(A, np.array([3., 6.0000001]))   # -> [1., 1.]
np.linalg.solve(A, np.array([3.0000001, 6.]))   # -> [7.0000001, -2.]

Три почти одинаковые правые части — три совершенно разных ответа. Изменение седьмого знака в правой части полностью меняет решение. При этом код отработал молча, вернул три красивых вектора, и ничто в выводе не говорит, что этим числам нельзя верить.

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

Особенно коварно, что «почти вырожденность» возникает сама собой из точно вырожденной задачи. Классический пример из ML: два признака — рост в сантиметрах и рост в дюймах. Математически второй столбец — это первый, делённый на $2{,}54$, матрица $X$ вырождена, и $X^TX$ обязана быть вырожденной. Но число $2{,}54$ не представимо точно в двоичной арифметике, и после умножения получается не точный ноль, а $\det(X^TX) \approx -1{,}1\cdot10^{-12}$:

X = np.array([[1., 2.54], [2., 5.08], [3., 7.62], [4., 10.16]])
G = X.T @ X
np.linalg.det(G)       # -1.08e-12, а не 0
np.linalg.cond(G)      # 1.07e+16

y = np.array([1., 2., 3., 4.])
np.linalg.solve(G, X.T @ y)      # -> [1.   , 0.   ]
np.linalg.inv(G) @ (X.T @ y)     # -> [0.125, 0.   ]

Исключения нет ни там, ни там. Два математически эквивалентных способа выдали разные ответы: $(1;\ 0)$ и $(0{,}125;\ 0)$. Оба неверны — точнее, оба являются произвольными точками из бесконечного множества решений вырожденной задачи, и какая именно точка получится, определяют ошибки округления. Если бы вместо 2.54 стояло 2.5, представимое в двоичной системе точно, был бы честный LinAlgError. То есть повезло бы больше.

Что с этим делать на практике. Проверять не определитель, а обусловленность: np.linalg.cond(A). Порог здоровья грубо такой — если $\operatorname{cond}(A) > 10^{12}$, ответу нельзя верить, и надо менять постановку задачи (регуляризация, отбор признаков, масштабирование), а не метод решения. И отдельно: если задача действительно вырождена — например, признаки линейно зависимы, — правильный ответ не «взять любое решение», а осознанно выбрать, какое из бесконечного множества решений тебе нужно; для этого существуют регуляризация и псевдообратная матрица, о которой речь пойдёт в следующих модулях курса.

Что до полностью вырожденного случая с нулевой правой частью, $Ax = 0$, — там всё интереснее: у такой системы всегда есть решение (нулевое), и вопрос в том, есть ли другие. Этому и посвящён следующий урок.

Решение систем и обусловленность в машинном обучении

Линейные системы в ML не где-то на периферии — они внутри почти каждого второго алгоритма. Ниже — карта того, где именно и какие инструменты для этого существуют.

Инструментарий: solve во всех библиотеках

NumPy. Базовый вызов — numpy.linalg.solve(A, b). Он принимает как вектор, так и матрицу правых частей: solve(A, B) решает все столбцы $B$ за одно разложение. Под капотом — LAPACK-функция dgesv: LU-разложение с частичным выбором ведущего элемента плюс две подстановки.

SciPy даёт то же самое, но с важной добавкой — флагом assume_a, которым ты сообщаешь решателю известные тебе свойства матрицы:

from scipy.linalg import solve

x = solve(A, b, assume_a='gen')   # общий случай (по умолчанию)
x = solve(A, b, assume_a='sym')   # матрица симметричная
x = solve(A, b, assume_a='pos')   # симметричная положительно определённая
x = solve(A, b, assume_a='her')   # эрмитова (для комплексных)

Флаг assume_a='pos' — самый ценный. Он говорит: «матрица симметрична, и все её диагональные миноры положительны». Для таких матриц существует разложение Холецкого — представление $A = LL^T$ через одну нижнюю треугольную матрицу вместо двух множителей LU. Оно вдвое дешевле обычного LU (примерно $\frac{n^3}{3}$ операций вместо $\frac{2n^3}{3}$), требует вдвое меньше памяти и не нуждается в выборе ведущего элемента — для положительно определённых матриц он не нужен, метод устойчив сам по себе. Технику разложения мы здесь не разбираем (это тема курса численных методов), но мотивация должна быть понятной: знание о структуре матрицы стоит денег, буквально — вдвое меньше времени за одну строчку в вызове.

Матрицы вида $X^TX$, ковариационные матрицы, матрицы Грама, гессианы выпуклых функций — все симметричны и положительно определены (или полуопределены). То есть почти в каждой ML-задаче, где возникает система, assume_a='pos' применим.

PyTorch. torch.linalg.solve(A, B) делает то же самое, но с двумя ключевыми отличиями. Первое: результат дифференцируем — градиент проходит сквозь решение системы, поэтому решатель можно поставить прямо внутрь вычислительного графа модели. Второе: он работает батчами.

Батчевое решение систем

Тензор формы (B, n, n) PyTorch (и NumPy тоже) интерпретирует как набор из $B$ независимых матриц $n\times n$, а не как одну трёхмерную матрицу. Вызов

A = torch.randn(64, 50, 50) + 50*torch.eye(50)   # 64 матрицы 50x50
b = torch.randn(64, 50, 1)                        # 64 правые части
X = torch.linalg.solve(A, b)                      # -> форма (64, 50, 1)

решает 64 разные системы за один вызов. Различие с solve(A, B) из предыдущего раздела принципиальное: там была одна матрица и много правых частей, здесь — много разных матриц, каждая со своей правой частью.

Зачем это нужно. В глубоком обучении данные почти всегда идут батчами: 64 объекта в мини-батче, у каждого своя маленькая линейная задача. Примеры: слой оптимизации внутри сети (differentiable optimization layers), где для каждого объекта решается своя квадратичная задача; нормализующие потоки, где обращается своя матрица преобразования; вероятностные модели, где у каждого объекта свой ковариационный блок. Наивный цикл for i in range(64): solve(A[i], b[i]) работал бы в десятки раз медленнее — не из-за арифметики, а из-за накладных расходов на 64 отдельных вызова и из-за невозможности загрузить GPU. Батчевый вызов передаёт всю пачку в LAPACK/cuBLAS одним куском, и ядра GPU перемалывают её параллельно.

Гауссовские процессы и jitter

Самый показательный сюжет про обусловленность во всём ML. Гауссовский процесс задаётся ковариационной функцией (ядром), и для предсказания нужно решить систему

где $K$ — матрица ковариаций между всеми обучающими точками, $K_{ij} = k(x_i, x_j)$.

Проблема в том, что матрица $K$ по своей природе плохо обусловлена. Если две обучающие точки близки, соответствующие строки $K$ почти одинаковы — а почти одинаковые строки, как мы уже видели, это прямая дорога к огромному $\operatorname{cond}$. И чем больше данных, тем гуще точки и тем хуже становится матрица.

Численный пример. Возьмём восемь точек на отрезке $[0;\ 0{,}7]$ с шагом $0{,}1$ и стандартное RBF-ядро $k(x, x') = \exp\!\left(-\frac{(x-x')^2}{2\ell^2}\right)$ с $\ell = 0{,}5$:

Десять знаков точности потеряны на ровном месте, притом что точек всего восемь. Наименьшее собственное число этой матрицы равно $2{,}3\cdot10^{-10}$ при наибольшем $6{,}70$ — их отношение и даёт обусловленность порядка $10^{10}$.

Стандартное лекарство — jitter (буквально «дрожание»): к диагонали прибавляют крошечное число $\varepsilon$:

Смотрим, что это даёт:

Добавка $10^{-6}$ — величина, физически бессмысленная на фоне значений ядра порядка единицы, — улучшила обусловленность в четыре тысячи раз. Механизм прозрачен: прибавление $\varepsilon$ к диагонали поднимает все собственные числа на $\varepsilon$, и наименьшее из них перестаёт быть «почти нулём». А поскольку $\operatorname{cond}$ в евклидовой норме для симметричной положительно определённой матрицы равен отношению наибольшего собственного числа к наименьшему, знаменатель перестаёт обнулять дробь.

Заметь, что это та же самая математика, что и в гребневой регрессии, где к $X^TX$ прибавляют $\lambda I$: одна операция, две интерпретации. В регрессии её называют регуляризацией и объясняют статистически (штраф за большие веса), в гауссовских процессах — jitter'ом и объясняют численно (шум наблюдений). Формула одна.

Фильтр Калмана

Фильтр Калмана — алгоритм, который в реальном времени оценивает скрытое состояние системы по зашумлённым измерениям: положение спутника по сигналам GPS, положение робота по показаниям одометрии и лидара, скрытая волатильность актива по котировкам. Он работает в навигации, в робототехнике, в трекинге объектов на видео.

На каждом шаге фильтра вычисляется так называемое усиление Калмана — матрица, определяющая, насколько сильно довериться новому измерению против предсказания модели:

где $S$ — ковариация невязки измерения. В учебниках это записывают через $S^{-1}$, а в реализациях, которые работают годами без сбоев, обращения нет: вместо $K = PH^TS^{-1}$ решают систему $S K^T = H P^T$ (матрица $S$ симметрична и положительно определена, так что подходит Холецкий). Причина — ровно наша: фильтр работает тысячи шагов подряд, ошибки накапливаются, и матрица $P$ со временем может потерять положительную определённость из-за округлений; явное обращение приближает этот момент.

Обусловленность матрицы признаков и масштабирование

И, наконец, самое повседневное применение всей этой теории — то, что ты делаешь в каждом пайплайне, часто не задумываясь.

Пусть в датасете два признака: площадь квартиры в квадратных метрах (значения порядка 50) и её цена в рублях (значения порядка 10 000 000). Столбцы матрицы $X$ различаются по масштабу в сотни тысяч раз. Что это делает с обусловленностью — на маленьком примере:

Теперь поделим второй столбец на 1000 — то есть просто сменим единицу измерения второго признака:

Обусловленность улучшилась в 357 раз от смены единиц измерения. Задача по сути та же, ответ пересчитывается тривиально, а численное поведение стало несопоставимо лучше.

Отсюда — прямое объяснение того, зачем нужен StandardScaler и его родня. Стандартизация признаков (вычесть среднее, поделить на стандартное отклонение) — это не косметика и не «требование алгоритма»; это операция, которая приводит столбцы к одному масштабу и тем самым снижает $\operatorname{cond}(X^TX)$. Следствия видны везде:

• линейные модели, решаемые через нормальные уравнения, дают более устойчивые коэффициенты (сам вывод нормальных уравнений и разбор мультиколлинеарности были в уроке про обратную матрицу, здесь мы смотрим на них только с численной стороны);

• градиентный спуск сходится быстрее: скорость сходимости для квадратичной задачи напрямую зависит от $\operatorname{cond}$ матрицы вторых производных, и при большом $\operatorname{cond}$ линии уровня функции потерь превращаются в вытянутый овраг, по которому градиент мечется от стенки к стенке вместо того, чтобы двигаться ко дну;

• регуляризация ($L_2$, она же ridge, она же weight decay) улучшает обусловленность по тому же механизму, что и jitter: прибавка $\lambda$ к диагонали поднимает наименьшее собственное число;

• методы второго порядка (Ньютон, L-BFGS, natural gradient) на каждом шаге решают систему с гессианом или его приближением — и их практическая пригодность целиком определяется обусловленностью этой матрицы.

Так что фраза «отмасштабируй признаки перед обучением» — это на самом деле «уменьши число обусловленности матрицы, с которой будет работать оптимизатор». Теперь у тебя есть инструмент, чтобы проверить это самому: посчитать np.linalg.cond(X.T @ X) до и после StandardScaler и посмотреть на разницу.

Практика: 30 заданий

Базовые задания (1–10)

Задание 1: Реши матричным методом систему $\begin{cases} 2x_1 + x_2 = 5 \\ x_1 + 3x_2 = 10 \end{cases}$ и проверь ответ подстановкой.

Задание 2: Реши матричным методом систему $\begin{cases} 3x_1 + 2x_2 = 7 \\ x_1 + 4x_2 = 9 \end{cases}$

Задание 3: Применим ли матричный метод к системе с матрицей $A = \begin{pmatrix} 1 & 2 & 3 \\ 2 & 4 & 6 \\ 1 & 0 & 1 \end{pmatrix}$? Если нет — объясни, что именно мешает.

Задание 4: Реши матричным методом систему $\begin{cases} 2x_1 - x_2 + x_3 = 3 \\ x_1 + 3x_2 - 2x_3 = 1 \\ 3x_1 + x_2 + x_3 = 8 \end{cases}$, если известно, что $A^{-1} = \dfrac{1}{9}\begin{pmatrix} 5 & 2 & -1 \\ -7 & -1 & 5 \\ -8 & -5 & 7 \end{pmatrix}$.

Задание 5: Вычисли все три нормы вектора $x = (1,\ -2,\ 2,\ 4)$.

Задание 6: Вычисли нормы $\|A\|_\infty$ и $\|A\|_1$ для матрицы $A = \begin{pmatrix} 1 & -5 & 2 \\ 3 & 0 & -1 \\ -2 & 4 & 4 \end{pmatrix}$.

Задание 7: Найди число обусловленности $\operatorname{cond}_\infty(A)$ для матрицы $A = \begin{pmatrix} 3 & 1 \\ 2 & 4 \end{pmatrix}$ и оцени, во сколько раз может усилиться относительная погрешность правой части.

Задание 8: Объясни, почему в системе $Ax = b$ (где $A$ размера $3\times3$) нельзя домножить обе части на $A^{-1}$ справа. Что конкретно не так с выражением $bA^{-1}$?

Задание 9: Реши матричным методом систему $\begin{cases} 4x_1 - x_2 = 9 \\ 2x_1 + 3x_2 = 13 \end{cases}$

Задание 10: Мастерская делает изделия двух типов. На одно изделие первого типа уходит 2 кг металла и 4 часа работы, на изделие второго типа — 3 кг металла и 1 час работы. За смену израсходовано 130 кг металла и 110 человеко-часов. Сколько изделий каждого типа сделали? Реши матричным методом.

Средние задания (11–20)

Задание 11: Реши матричным методом систему $\begin{cases} 3x_1 + x_2 + 2x_3 = 3 \\ x_1 + 2x_2 + x_3 = 4 \\ 2x_1 + x_2 + 3x_3 = 1 \end{cases}$

Задание 12: Реши матричное уравнение $AX = B$ с двумя правыми частями: $A = \begin{pmatrix} 5 & 2 \\ 3 & 4 \end{pmatrix}$, $B = \begin{pmatrix} 16 & 1 \\ 18 & -5 \end{pmatrix}$.

Задание 13: Реши $AX = B$ с тремя правыми частями: $A = \begin{pmatrix} 1 & 0 & 2 \\ 2 & -1 & 3 \\ 4 & 1 & 8 \end{pmatrix}$, $B = \begin{pmatrix} 1 & 2 & 4 \\ 0 & 0 & 8 \\ 6 & 11 & 15 \end{pmatrix}$, если $A^{-1} = \begin{pmatrix} -11 & 2 & 2 \\ -4 & 0 & 1 \\ 6 & -1 & -1 \end{pmatrix}$.

Задание 14: При каких значениях параметра $a$ к системе $\begin{cases} x_1 + 2x_2 = 3 \\ 3x_1 + a\,x_2 = 5 \end{cases}$ применим матричный метод? Найди решение при $a = 4$.

Задание 15: При каких значениях параметра $a$ применим матричный метод к системе с матрицей $A = \begin{pmatrix} 1 & 1 & 1 \\ 1 & 2 & a \\ 1 & 4 & a^2 \end{pmatrix}$?

Задание 16: Найди $\operatorname{cond}_\infty(A)$ для матрицы $A = \begin{pmatrix} 1 & 2 \\ 2 & 4{,}001 \end{pmatrix}$ и оцени, сколько верных десятичных знаков останется в решении, если правая часть известна с относительной погрешностью $10^{-6}$.

Задание 17: Известно, что $\operatorname{cond}(A) = 250$, а правая часть измерена с относительной погрешностью 0,2 %. Какова гарантированная оценка относительной погрешности решения? Можно ли доверять первому знаку после запятой в ответе, если решение имеет порядок единицы?

Задание 18: Сравни число арифметических операций для системы $6\times6$ тремя методами: Гауссом ($\frac{2n^3}{3} + 2n^2$), матричным ($2n^3 + 2n^2$) и правилом Крамера с разложением определителей по строке ($\approx (n+1)\cdot n!$).

Задание 19: Фабрика выпускает три вида продукции, расходуя три ресурса; матрица расхода $A = \begin{pmatrix} 1 & 2 & 1 \\ 2 & 1 & 3 \\ 1 & 1 & 1 \end{pmatrix}$ (строка — ресурс, столбец — продукт). Известно, что $A^{-1} = \begin{pmatrix} -2 & -1 & 5 \\ 1 & 0 & -1 \\ 1 & 1 & -3 \end{pmatrix}$. Найди планы выпуска для двух сценариев запасов: $b^{(1)} = (80,\ 130,\ 60)$ и $b^{(2)} = (80,\ 170,\ 70)$.

Задание 20: Реши матричным методом систему $4\times4$: $\begin{cases} x_1 + x_2 + x_3 + x_4 = 10 \\ x_2 + x_3 + x_4 = 9 \\ x_3 + x_4 = 7 \\ x_4 = 4 \end{cases}$

Продвинутые задания (21–30)

Задание 21: Для матрицы Гильберта $H = \begin{pmatrix} 1 & \frac12 & \frac13 \\[3pt] \frac12 & \frac13 & \frac14 \\[3pt] \frac13 & \frac14 & \frac15 \end{pmatrix}$ известно, что $H^{-1} = \begin{pmatrix} 9 & -36 & 30 \\ -36 & 192 & -180 \\ 30 & -180 & 180 \end{pmatrix}$. Найди $\operatorname{cond}_\infty(H)$ и оцени, какой будет относительная погрешность решения, если правую часть округлить до трёх знаков после запятой (решение системы — вектор $(1,1,1)$).

Задание 22: Дана система $\begin{cases} x_1 + x_2 = 2 \\ x_1 + 1{,}001\,x_2 = 2{,}001 \end{cases}$. Найди решение. Затем измени второй свободный член на $2{,}002$ и найди новое решение. Посчитай $\operatorname{cond}_\infty(A)$ и фактическое усиление относительной погрешности (в норме $\|\cdot\|_\infty$).

Задание 23: Реши буквенную систему $\begin{cases} k\,x_1 + x_2 = 1 \\ x_1 + k\,x_2 = 1 \end{cases}$ матричным методом. При каких $k$ метод применим? Что происходит с решением при $k \to -1$?

Задание 24: Докажи, что $\operatorname{cond}(A) \ge 1$ для любой невырожденной матрицы, и что $\operatorname{cond}(\alpha A) = \operatorname{cond}(A)$ для любого числа $\alpha \ne 0$. Какие матрицы имеют $\operatorname{cond} = 1$?

Задание 25: Что выведет этот код и почему?

import numpy as np
A = np.array([[2., 4.], [3., 6.]])
b = np.array([10., 15.])
x = np.linalg.solve(A, b)
print(x)

numpy.linalg.LinAlgError: Singular matrix

Задание 26: А что выведет этот код? Объясни результат.

import numpy as np
A = np.array([[1., 2.], [2., 4.0000001]])
print(np.linalg.cond(A))
print(np.linalg.solve(A, np.array([3., 6.])))
print(np.linalg.solve(A, np.array([3., 6.0000001])))

250000000.80668727
[ 3. -0.]
[1. 1.]

Задание 27: Вычисления идут в типе float64 (около 16 верных десятичных знаков). Сколько верных знаков останется в решении системы при $\operatorname{cond}(A) = 10^{7}$? А какое максимальное $\operatorname{cond}$ допустимо, если тебе нужны хотя бы 3 верных знака в ответе?

Задание 28: Реши матричное уравнение $AX = E$ для $A = \begin{pmatrix} 2 & 1 \\ 5 & 3 \end{pmatrix}$, где $E$ — единичная матрица $2\times2$. Что представляет собой найденная матрица $X$? Как этот результат связывает решение систем и обращение матриц?

Задание 29: Признаки в датасете измерены в разных единицах, и матрица получилась такой: $A = \begin{pmatrix} 1 & 1000 \\ 2 & 1500 \end{pmatrix}$. Посчитай $\operatorname{cond}_\infty(A)$. Затем измени единицу измерения второго признака, поделив второй столбец на 1000, и посчитай обусловленность заново. Во сколько раз стало лучше?

Задание 30 (капстоун): Дана система

Выполни полный разбор: (а) проверь применимость матричного метода; (б) найди $A^{-1}$; (в) реши систему и проверь подстановкой; (г) посчитай $\operatorname{cond}_\infty(A)$; (д) оцени погрешность решения, если правая часть измерена с точностью 1 %; (е) объясни, как ты решал бы эту систему в коде и почему.

from scipy.linalg import solve
x = solve(A, b, assume_a='pos')   # разложение Холецкого

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

Ошибка 1. Записывают решение как $x = bA^{-1}$ или пытаются «поделить $b$ на $A$».

Как выглядит: студент по аналогии с числами пишет $x = b / A$ или ставит $A^{-1}$ справа от правой части. Почему возникает: перенос школьной привычки, где $ax = b \Rightarrow x = b/a$, и порядок множителей не важен. Как правильно: матричное умножение некоммутативно, а деления на матрицу вообще не существует. Единственная корректная запись — $x = A^{-1}b$, с $A^{-1}$ слева. В системе $Ax = b$ выражение $bA^{-1}$ даже не определено по размерам: $(n\times1)\cdot(n\times n)$ — такого произведения не бывает.

Ошибка 2. Применяют матричный метод, не проверив определитель.

Как выглядит: человек сразу бросается считать союзную матрицу, тратит десять минут, доходит до деления на $\det A$ и обнаруживает там ноль. Почему возникает: проверка кажется формальностью «для галочки». Как правильно: $\det A \ne 0$ — это критерий применимости, а не ритуал. Считать определитель нужно первым действием, до всех остальных вычислений. Если он равен нулю, метод не применим вовсе, и нужно переходить к методу Гаусса.

Ошибка 3. Считают, что $\det A = 0$ означает «у системы нет решений».

Как выглядит: получив нулевой определитель, пишут в ответе «система несовместна». Почему возникает: смешение двух разных утверждений — «нет единственного решения» и «нет решений». Как правильно: $\det A = 0$ означает только то, что единственного решения нет. Дальше возможны оба варианта: решений может не быть совсем, а может быть бесконечно много. Различить их умеет метод Гаусса (или сравнение рангов по теореме Кронекера — Капелли), а матричный метод — нет.

Ошибка 4. Путают маленький определитель с плохой обусловленностью.

Как выглядит: «определитель $10^{-6}$ — значит, задача плохо обусловлена». Почему возникает: интуиция «маленький определитель — близко к вырожденности» верна только качественно. Как правильно: определитель зависит от масштаба матрицы, а обусловленность — нет. У матрицы $0{,}01\cdot E$ размера $3\times3$ определитель равен $10^{-6}$, а $\operatorname{cond} = 1$ — она идеальна. И наоборот, у матрицы с определителем порядка единицы обусловленность может быть $10^{10}$. Индикатор численных проблем — это $\operatorname{cond}(A)$, а не $\det A$.

Ошибка 5. Считают число обусловленности точным предсказанием ошибки.

Как выглядит: «$\operatorname{cond} = 1000$, значит ошибка будет ровно в тысячу раз больше». Почему возникает: знак $\le$ в формуле проскакивают взглядом. Как правильно: $\operatorname{cond}$ даёт верхнюю границу для наихудшего направления возмущения. Конкретная ошибка может оказаться сильно меньше — мы видели это на матрице Гильберта, где возмущение вдоль самого вектора $b$ дало усиление ровно 1 вместо возможных 15 514. Обусловленность отвечает на вопрос «насколько плохо может быть», а не «насколько плохо будет».

Ошибка 6. В коде пишут np.linalg.inv(A) @ b.

Как выглядит: перевод формулы $x = A^{-1}b$ в код буква в букву. Почему возникает: формула из учебника выглядит как инструкция к действию. Как правильно: np.linalg.solve(A, b). Учебная формула описывает что найти, а не как считать. Явная обратная матрица требует втрое больше операций, даёт ошибку в 15–70 раз больше и невязку на восемь порядков хуже — при абсолютно том же математическом результате.

Ошибка 7. Проверяют вырожденность через if np.linalg.det(A) == 0.

Как выглядит: защитная проверка перед решением системы. Почему возникает: прямой перевод математического критерия в код. Как правильно: в арифметике с плавающей точкой определитель почти никогда не равен ровно нулю — вспомни пример с ростом в сантиметрах и дюймах, где $\det = -1{,}08\cdot10^{-12}$ вместо нуля. Сравнивать вещественное число с нулём через == бессмысленно. Проверять надо обусловленность: np.linalg.cond(A) > 1e12 — и это тоже эвристика, а не строгий критерий.

Ошибка 8. Забывают вынести $\frac{1}{\det A}$ за скобку при ручном счёте.

Как выглядит: каждый элемент союзной матрицы делится на определитель, дальше идёт умножение на $b$ в дробях, и на третьем шаге теряется знаменатель. Почему возникает: буквальное следование формуле $A^{-1} = \frac{1}{\det A}\operatorname{adj}(A)$. Как правильно: держать $\frac{1}{\det A}$ снаружи, всё умножение делать в целых числах и разделить один раз в самом конце. Это не только быстрее, но и сильно снижает шанс арифметической ошибки — сравни примеры 1 и 2 из теоретической части.

Ошибка 9. Строят $A^{-1}$ ради одной-единственной системы.

Как выглядит: «мне надо решить $Ax = b$, значит сначала найду обратную матрицу». Почему возникает: формула $x = A^{-1}b$ читается как последовательность шагов. Как правильно: для одной правой части это втрое больше работы, чем нужно. Обратная матрица окупается начиная примерно с трёх правых частей — и даже тогда честнее один раз разложить матрицу (LU или Холецкий) и переиспользовать разложение.

Ошибка 10. Думают, что смена алгоритма спасёт плохо обусловленную задачу.

Как выглядит: «solve дал плохой ответ, попробую другой решатель / другую библиотеку / другой язык». Почему возникает: смешение проблемы метода и проблемы задачи. Как правильно: если $\operatorname{cond}(A) = 10^{14}$, то ни один алгоритм не даст больше двух верных знаков — информации нет в самих данных. Менять надо постановку: регуляризовать, отбросить почти зависимые признаки, перемасштабировать переменные, перейти к более высокой точности арифметики. Смена решателя лечит только вторую часть проблемы, ту, что от алгоритма.

Главное запомнить

• Матричный метод: если $A$ квадратная и $\det A \ne 0$, то система $Ax = b$ имеет единственное решение $x = A^{-1}b$.

• Домножать на $A^{-1}$ нужно слева, потому что умножение матриц некоммутативно; справа выражение $bA^{-1}$ вообще не определено по размерам. Общее правило: множитель убирают с той стороны, с которой он стоит.

• Проверка $\det A \ne 0$ — первое действие, а не формальность: при нулевом определителе метод не применим (что не означает отсутствия решений — их может быть бесконечно много).

• Несколько систем с одной матрицей и разными правыми частями записываются как $AX = B$ и решаются одним умножением $X = A^{-1}B$; столбцы работают независимо друг от друга.

• Явная обратная матрица окупается примерно с трёх правых частей; для одной системы она стоит $2n^3$ операций против $\frac{2n^3}{3}$ у метода Гаусса — втрое дороже при том же результате.

• Норма вектора — его «длина» ($\|x\|_2$, $\|x\|_1$, $\|x\|_\infty$); норма матрицы — максимальный коэффициент растяжения, $\|A\| = \max_{x\ne0}\frac{\|Ax\|}{\|x\|}$. Руками удобно считать $\|A\|_\infty$ как максимальную сумму модулей по строке.

• Число обусловленности $\operatorname{cond}(A) = \|A\|\cdot\|A^{-1}\|$ показывает, во сколько раз может усилиться относительная погрешность: $\dfrac{\|\delta x\|}{\|x\|} \le \operatorname{cond}(A)\cdot\dfrac{\|\delta b\|}{\|b\|}$.

• Всегда $\operatorname{cond}(A) \ge 1$, и $\operatorname{cond}$ не меняется при умножении матрицы на число. Индикатор численных проблем — это $\operatorname{cond}$, а не маленький определитель.

• Правило потери знаков: теряется примерно $\lg\operatorname{cond}(A)$ верных десятичных знаков. В float64 (16 знаков) при $\operatorname{cond} = 10^{16}$ не остаётся ни одного.

• В коде: numpy.linalg.solve(A, b), а не numpy.linalg.inv(A) @ b. Причины — точность (ошибка в десятки раз меньше, невязка на порядки лучше), скорость (втрое) и память (обратная к разреженной матрице плотная).

• Вырожденная матрица падает с LinAlgError: Singular matrix — и это хороший исход. Почти вырожденная не падает, а молча возвращает мусор — это исход опасный.

• Плохую обусловленность лечат не сменой алгоритма, а сменой постановки: масштабированием признаков, регуляризацией (добавка $\varepsilon I$ или $\lambda I$ к диагонали), удалением почти зависимых столбцов.

Связь с другими темами курса

Что нужно было знать до этого урока

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

Из урока про обратную матрицу (161) пришёл сам объект $A^{-1}$: критерий существования $\det A \ne 0$, способы построения через союзную матрицу и методом Гаусса, а также правило работы с матричными уравнениями $AX = B$ и $XA = B$. В этом уроке мы не искали обратную матрицу заново — мы применяли её как готовый инструмент и смотрели на неё с новой стороны: как на оператор, переводящий правую часть в решение.

Из уроков про определители (158–160) пришло умение быстро проверять невырожденность и считать союзную матрицу через алгебраические дополнения — без этого не сделать ни одного шага матричного метода вручную.

Из урока про ранг матрицы (162) — понимание того, что вырожденность есть падение ранга ниже $n$, то есть наличие «лишних», зависимых уравнений.

Из уроков про системы линейных уравнений (163), метод Гаусса (164) и правило Крамера (165) — вся рамка задачи: матричная запись $Ax = b$, теорема Кронекера — Капелли, критерии совместности и определённости, алгоритм Гаусса с его прямым и обратным ходом и оценкой $O(n^3)$, крамеровские формулы и их стоимость. Именно на этом фоне матричный метод занимает своё место третьим — и мы смогли сравнить все три честно, с числами.

Что изучить дальше

Урок 167 «Однородные системы» — прямое продолжение. Матричный метод работает только при $\det A \ne 0$, а следующий урок разбирает ровно противоположный случай: систему $Ax = 0$ с вырожденной матрицей, где решений бесконечно много, и их множество надо описать целиком через фундаментальную систему решений.

Дальше в курсе линейной алгебры понятия этого урока раскроются глубже. Векторные пространства (168–170) дадут строгий язык для нормы: нормированное пространство, метрика, полнота. Собственные числа и векторы (172) объяснят, откуда берётся обусловленность симметричных матриц: для них $\operatorname{cond}_2$ — это в точности отношение наибольшего собственного числа к наименьшему по модулю. Квадратичные формы (174) и евклидово пространство (175) дадут точное определение положительной определённости — того самого свойства, которое включает флаг assume_a='pos' и разложение Холецкого.

За пределами этого курса тебя ждут LU- и QR-разложения, сингулярное разложение (SVD) и псевдообратная матрица Мура — Пенроуза — инструменты, которые аккуратно решают всё то, где матричный метод пасует: прямоугольные системы, вырожденные матрицы, задачи наименьших квадратов. Именно через сингулярные числа даётся самое употребительное определение обусловленности, $\operatorname{cond}_2(A) = \sigma_{\max}/\sigma_{\min}$.

Где это нужно в жизни

💻 Программирование. Любая работа с численными библиотеками требует знать разницу между solve и inv. В компьютерной графике решают системы для обращения матриц трансформации камеры и для инверсной кинематики скелетной анимации. В физических движках игр каждый кадр решается система для расчёта импульсов в контактах. В решателях уравнений в частных производных (гидродинамика, теплопроводность, прочность) системы бывают размером в миллионы неизвестных — и там всё держится на разреженных решателях, которые никогда не строят обратную матрицу.

🤖 ML/AI. Гауссовские процессы (решение системы с ковариационной матрицей и jitter), гребневая регрессия, методы оптимизации второго порядка (Ньютон, L-BFGS, natural gradient), нормализующие потоки, дифференцируемые слои оптимизации, фильтр Калмана в трекинге объектов. Понимание обусловленности напрямую объясняет, зачем нужны нормализация признаков, batch normalization, weight decay и градиентный клиппинг.

📊 Data Science. Диагностика мультиколлинеарности через $\operatorname{cond}(X^TX)$ и VIF, объяснение нестабильных коэффициентов регрессии, выбор параметра регуляризации, задачи балансировки и распределения ресурсов, где матрица затрат одна, а сценариев много.

🔬 Наука. В геодезии и астрономии — уравнивание сетей измерений методом наименьших квадратов (та самая задача, ради которой Гаусс развивал свой метод). В химии — балансировка уравнений реакций и расчёт равновесных концентраций. В структурной механике — расчёт напряжений в фермах и конструкциях. В томографии — реконструкция изображения по проекциям, задача классически плохо обусловленная, где без регуляризации не обойтись.

💰 Финансы. Модель «затраты — выпуск» Леонтьева, где матрица $(E - A)^{-1}$ вычисляется явно ради своих элементов. Оптимизация портфеля по Марковицу, где решается система с ковариационной матрицей доходностей, — и она почти всегда плохо обусловлена, потому что доходности активов сильно коррелированы, поэтому её обязательно регуляризуют (shrinkage-оценки Ледуа — Вольфа — это буквально добавка к диагонали). Калибровка моделей ценообразования опционов, актуарные расчёты.

Интересные факты

• Термин «число обусловленности» ввёл Алан Тьюринг в статье 1948 года «Rounding-off Errors in Matrix Processes» — той же самой, где он ввёл и понятие LU-разложения (он называл его «triangular resolution»). Тьюринг занимался этим не из теоретического интереса, а потому что проектировал реальную вычислительную машину ACE и ему нужно было понимать, каким результатам можно верить.

• Матрицу Гильберта, ставшую эталонным примером плохой обусловленности, Давид Гильберт нашёл в 1894 году — за полвека до появления электронных компьютеров и задолго до того, как слово «обусловленность» вообще возникло. Её число обусловленности растёт примерно как $e^{3{,}5n}$: для $n = 6$ это $10^{7}$, для $n = 10$ — $10^{13}$, для $n = 13$ — уже $10^{18}$, то есть в стандартной машинной арифметике решение системы с ней теряет всякий смысл.

• В документации MATLAB функция inv сопровождается прямым предупреждением не использовать её для решения систем, а рекомендованный оператор \ (backslash) выбирает алгоритм автоматически: он проверяет, треугольная ли матрица, симметричная ли, положительно определённая ли, разреженная ли, — и подбирает подходящий метод из десятка вариантов. Одна из причин, по которой MATLAB десятилетиями оставался стандартом инженерных расчётов, — именно этот оператор.

• Джеймс Уилкинсон, автор обратного анализа ошибок, во время Второй мировой войны работал в баллистической лаборатории, а затем с Тьюрингом над проектом ACE. Его знаменитая фраза о том, что за годы работы он не встретил ни одной практической задачи, где метод Гаусса с частичным выбором ведущего элемента оказался бы численно неустойчивым, до сих пор остаётся ориентиром: теоретическая оценка роста элементов у этого метода экспоненциальная, а на практике он ведёт себя превосходно почти всегда. Разрыв между худшим случаем и практикой здесь один из самых больших во всей вычислительной математике.

• Матрица полных затрат $(E - A)^{-1}$ из модели Василия Леонтьева — один из редких случаев, когда обратную матрицу вычисляют именно как объект, а не как способ решить систему: её элементы читаются экономически, как полные (прямые плюс косвенные) затраты одной отрасли на единицу продукции другой. Леонтьев получил за эту модель Нобелевскую премию по экономике в 1973 году, а первые расчёты по экономике США он проводил в 1940-х на матрицах порядка 40 — обращая их вручную и на электромеханических табуляторах, что заняло месяцы работы.

Лайфхаки и полезные трюки

1. Определитель первым действием, а множитель $\frac{1}{\det A}$ — за скобку до конца. Тридцать секунд на $\det A$ экономят десять минут вычислений, которые при $\det A = 0$ обязаны провалиться; побочная польза — если определитель оказался равен $\pm1$, обратная матрица будет целочисленной, и дальше можно считать вообще без дробей. А сам множитель $\frac{1}{\det A}$ не размазывай по элементам союзной матрицы: перемножь всё в целых числах и раздели один раз в финале. Пример: вместо $\begin{pmatrix} -4/5 & -3/5 & 1 \\ 7/5 & 4/5 & -1 \\ 1 & 1 & -1\end{pmatrix}\begin{pmatrix} 5 \\ 0 \\ 6\end{pmatrix}$ считай $\frac15\begin{pmatrix} -4 & -3 & 5 \\ 7 & 4 & -5 \\ 5 & 5 & -5\end{pmatrix}\begin{pmatrix} 5 \\ 0 \\ 6\end{pmatrix} = \frac15\begin{pmatrix} 10 \\ 5 \\ -5\end{pmatrix}$.

2. Проверяй ответ подстановкой в исходную систему, а не в $A^{-1}$. Подстановка в исходные уравнения ловит ошибку на любом этапе — и в определителе, и в союзной матрице, и в финальном умножении. Проверка $A A^{-1} = E$ ловит ошибку только в обратной матрице. При этом подстановка занимает секунд двадцать: три скалярных произведения.

3. Для симметричной матрицы считай только шесть чисел вместо девяти. Союзная матрица симметричной матрицы тоже симметрична, значит $A_{12} = A_{21}$, $A_{13} = A_{31}$, $A_{23} = A_{32}$. Для $3\times3$ достаточно посчитать три диагональных дополнения и три наддиагональных. Симметричные матрицы в задачах встречаются постоянно — ковариационные, Грама, $X^TX$.

4. Оценивай обусловленность на глаз по «почти пропорциональности» строк. Если две строки матрицы отличаются в третьем-четвёртом знаке — жди $\operatorname{cond}$ порядка тысяч и выше. Геометрический образ: почти параллельные прямые (или почти совпадающие плоскости), точка пересечения которых мечется от малейшего сдвига. Наоборот, диагональное преобладание (когда $|a_{ii}|$ больше суммы модулей остальных элементов строки) — надёжный признак хорошей обусловленности.

5. В коде — три строки диагностики перед решением любой серьёзной системы.

print(A.shape)                  # квадратная ли
print(np.linalg.cond(A))        # можно ли верить ответу
x = np.linalg.solve(A, b)
print(np.linalg.norm(A @ x - b) / np.linalg.norm(b))   # невязка

Обусловленность выше $10^{12}$ — сигнал менять постановку задачи. Невязка выше $10^{-12}$ при разумном $\operatorname{cond}$ — сигнал, что что-то не так с самим кодом.

6. Помни правило потери знаков как устный калькулятор. $\lg\operatorname{cond}(A)$ — сколько десятичных знаков съедено. float64 даёт около 16 знаков, float32 — около 7. Отсюда важное следствие для ML: если ты обучаешь модель в float32 (а на GPU это норма) и где-то решаешь систему с $\operatorname{cond} = 10^{6}$, у тебя останется в лучшем случае один верный знак — в такой ситуации либо переходи на float64 для этого конкретного места, либо регуляризуй матрицу.

7. Добавка $\varepsilon I$ к диагонали — универсальное первое средство. Одна и та же операция называется по-разному в разных областях: jitter в гауссовских процессах, ridge/$L_2$/weight decay в регрессии и нейросетях, shrinkage в оценке ковариаций, демпфирование Левенберга — Марквардта в оптимизации. Механизм везде один: поднять наименьшее собственное число матрицы и тем самым сбить $\operatorname{cond}$. Начинать разумно с $\varepsilon$ порядка $10^{-6}$ от масштаба диагонали и увеличивать, пока обусловленность не станет приемлемой.

Матричный метод — редкий случай темы, где главный урок оказывается не в формуле, а в оговорках к ней. Формула $x = A^{-1}b$ проста, честна и абсолютно верна математически; всё интересное начинается там, где к ней добавляются вопросы «а какой ценой», «а с какой точностью» и «а что будет, если данные чуть-чуть другие». Именно эти вопросы отличают человека, который умеет решать системы, от человека, который умеет решать реальные задачи с системами внутри. И число обусловленности — тот единственный инструмент, который на все три вопроса отвечает одним числом.

В следующем уроке мы вернёмся к чистой алгебре и разберём случай, который матричный метод честно обходит стороной: систему $Ax = 0$ с вырожденной матрицей. У неё решений либо одно (тривиальное), либо сразу бесконечно много, и описать их все — задача не менее интересная, чем найти единственное.