Метод Гаусса: движок за каждым численным решением

Метод Гаусса работает внутри каждой библиотеки линейной алгебры - numpy, LAPACK, scipy. Когда Python вызывает `np.linalg.solve`, под капотом происходит именно гауссово исключение с частичным выбором ведущего элемента. Это алгоритм, которому 200 лет и который запущен на миллиардах устройств прямо сейчас.

• numpy.linalg.solve: внутри - LAPACK dgesv, который реализует метод Гаусса с LU
• Компьютерная графика: обратная матрица трансформации - через метод Гаусса
• Физика: системы уравнений МКЭ (конечных элементов) - миллионы переменных
• Экономика: модель общего равновесия - большая разреженная система
• Криптография: решение систем над полем GF(2) - атаки на коды с исправлением ошибок

Метод Гаусса: движок за каждым численным решением

**`np.linalg.solve(A, b)` - одна строка кода, которая решает любую систему линейных уравнений.** Внутри неё - LU-разложение, которое является прямым потомком алгоритма, описанного Гауссом в 1809 году. Когда scikit-learn обучает линейную регрессию на нормальном уравнении, когда OpenCV находит гомографию камеры, когда физический симулятор считает силы в системе тел - везде работает один и тот же механизм: последовательное исключение переменных. Метод Гаусса - это не "учебный приём", это производственный движок.

**О чём этот урок на самом деле**: не о том, как вручную решать системы 3×3. А о том, как алгоритм превращает запутанную систему в треугольную форму, откуда ответ читается напрямую - и почему именно этот алгоритм живёт внутри каждой научной библиотеки от NumPy до MATLAB.

Задача: найти x, y, z одновременно

Задача: найти x, y, z одновременно

Система трёх уравнений с тремя неизвестными - классический случай. Запишем её в матричном виде A·x = b:

СИСТЕМА: 2x + 1y - 1z = 8 4x - 6y + 0z = -2 -2x + 7y + 2z = 1 МАТРИЧНЫЙ ВИД A·x = b: A = [ 2 1 -1 ] x = [x] b = [ 8] [ 4 -6 0 ] [y] [-2] [-2 7 2 ] [z] [ 1] ЗАДАЧА: найти вектор x = [x, y, z] такой, что A·x = b ЭТО ОДНА И ТА ЖЕ ЗАДАЧА: - "решить систему уравнений" - "найти x = A⁻¹·b" (когда A обратима) - "подобрать коэффициенты под данные" (регрессия)

Расширенная матрица: система в одном объекте

Расширенная матрица: система в одном объекте

Первый шаг - записать систему как **расширенную матрицу** [A | b]. Правый столбец - правые части уравнений. Все дальнейшие операции - над строками этой матрицы:

[ 2 1 -1 | 8 ] [ 4 -6 0 | -2 ] [-2 7 2 | 1 ] КАЖДАЯ СТРОКА = одно уравнение КАЖДЫЙ СТОЛБЕЦ = один коэффициент ЦЕЛЬ: привести к такому виду: [ 1 * * | * ] <- треугольная (ступенчатая) форма [ 0 1 * | * ] [ 0 0 1 | * ] Из неё x, y, z читаются снизу вверх за 3 подстановки.

Три элементарных операции (которые не меняют решение)

Три элементарных операции (которые не меняют решение)

Над строками расширенной матрицы разрешены три операции. Каждая из них эквивалентна алгебраическому преобразованию системы - решение при этом не меняется:

**Pivot** - ведущий элемент строки. На каждом шаге прямого хода он обнуляет всё под собой. Выбор правильного pivot (максимального по модулю) - ключ к численной устойчивости. Именно это делает `partial pivoting` в numpy.

Прямой ход: шаг за шагом

Прямой ход: шаг за шагом

Прямой ход - это последовательное обнуление элементов **ниже главной диагонали**. Разберём по шагам на нашей системе:

ИСХОДНАЯ МАТРИЦА: [ 2 1 -1 | 8 ] [ 4 -6 0 | -2 ] [-2 7 2 | 1 ] ШАГ 1: pivot = 2 (строка 1, столбец 1) R2 ← R2 - (4/2)·R1 = R2 - 2·R1 R3 ← R3 - (-2/2)·R1 = R3 + 1·R1 [ 2 1 -1 | 8 ] [ 0 -8 2 | -18] <- 4-2·2=0, -6-2·1=-8, 0-2·(-1)=2, -2-2·8=-18 [ 0 8 1 | 9 ] <- -2+2=0, 7+1=8, 2+(-1)=1, 1+8=9 ШАГ 2: pivot = -8 (строка 2, столбец 2) R3 ← R3 - (8/-8)·R2 = R3 + 1·R2 [ 2 1 -1 | 8 ] [ 0 -8 2 | -18] [ 0 0 3 | -9 ] <- 8+(-8)=0, 1+2=3, 9+(-18)=-9 РЕЗУЛЬТАТ: треугольная форма - готова к обратному ходу

Обратный ход: читаем ответ снизу вверх

Обратный ход: читаем ответ снизу вверх

Треугольная форма позволяет вычислить переменные поочерёдно - начиная с последней строки, где осталась только одна неизвестная:

ТРЕУГОЛЬНАЯ СИСТЕМА: 2x + y - z = 8 - 8y + 2z = -18 3z = -9 СТРОКА 3: 3z = -9 => z = -3 СТРОКА 2: -8y + 2·(-3) = -18 -8y = -18 + 6 = -12 y = 3/2 СТРОКА 1: 2x + (3/2) - (-3) = 8 2x = 8 - 3/2 - 3 = 7/2 x = 7/4 ОТВЕТ: x = 1.75, y = 1.5, z = -3 ПРОВЕРКА: 2·(1.75) + 1·(1.5) - 1·(-3) = 3.5 + 1.5 + 3 = 8 ✓ 4·(1.75) - 6·(1.5) + 0·(-3) = 7 - 9 = -2 ✓ -2·(1.75) + 7·(1.5) + 2·(-3) = -3.5 + 10.5 - 6 = 1 ✓

Метод Гаусса в numpy: как это выглядит в коде

Метод Гаусса в numpy: как это выглядит в коде

**LU = метод Гаусса в матричном виде**: L хранит множители, которые использовались при прямом ходе. U - это результат прямого хода. Разложение A = P·L·U делается один раз, а потом система решается для любого вектора b за O(n²) вместо O(n³). Именно поэтому numpy сначала строит LU, а потом решает - даже для одной системы.

ML-применение №1: нормальное уравнение линейной регрессии

ML-применение №1: нормальное уравнение линейной регрессии

Аналитическое решение линейной регрессии - **нормальное уравнение** - это система линейных уравнений, которую решают методом Гаусса (через LU или через Cholesky для симметричного XᵀX):

ЗАДАЧА РЕГРЕССИИ: Минимизировать ||Xθ - y||² АНАЛИТИЧЕСКОЕ РЕШЕНИЕ: θ = (XᵀX)⁻¹ · Xᵀy ЭТО ЭКВИВАЛЕНТНО РЕШЕНИЮ СИСТЕМЫ: (XᵀX) · θ = Xᵀy A · θ = b где A = XᵀX, b = Xᵀy METOD ГАУССА (через LU-разложение) даёт θ без явного вычисления обратной матрицы.

ML-применение №2: метод Гаусса-Жордана для обратной матрицы

ML-применение №2: метод Гаусса-Жордана для обратной матрицы

Расширение метода Гаусса - **Гаусс-Жордан** - доводит матрицу до единичной не только снизу диагонали, но и сверху. Если приписать справа единичную матрицу, после редукции правая часть станет обратной матрицей. Именно этот алгоритм лежит в основе `np.linalg.inv`:

НАЙТИ ОБРАТНУЮ ДЛЯ A = [2 1]: [1 1] ДОПИСАТЬ ЕДИНИЧНУЮ МАТРИЦУ: [ 2 1 | 1 0 ] [ 1 1 | 0 1 ] R1 ← R1 / 2: [ 1 0.5 | 0.5 0 ] [ 1 1 | 0 1 ] R2 ← R2 - R1: [ 1 0.5 | 0.5 0 ] [ 0 0.5 | -0.5 1 ] R2 ← R2 / 0.5: [ 1 0.5 | 0.5 0 ] [ 0 1 | -1 2 ] R1 ← R1 - 0.5·R2: [ 1 0 | 1 -1 ] [ 0 1 | -1 2 ] A⁻¹ = [ 1 -1] ПРОВЕРКА: [2 1] · [ 1 -1] = [1 0] = I ✓ [-1 2] [1 1] [-1 2] [0 1]

**Метод Гаусса-Жордана находит обратную матрицу** - и именно он же используется при обучении нейросетей в редких случаях, когда нужна явная обратная (например, в second-order optimizers типа K-FAC). Обычно это заменяется итеративными методами, но понимание Гаусса-Жордана - основа.

Численная устойчивость и partial pivoting

Численная устойчивость и partial pivoting

Наивный метод Гаусса ломается, если на диагонали оказывается нуль или очень маленькое число: деление на него вносит гигантскую погрешность. **Partial pivoting** - перестановка строк перед каждым шагом, чтобы ведущим элементом стал максимальный в столбце:

БЕЗ PIVOTING (наивно): Система: 0.001·x + 1·y = 1 1·x + 1·y = 2 Pivot = 0.001 -> множитель = 1 / 0.001 = 1000 R2 ← R2 - 1000·R1 -> числа растут, накапливается ошибка С PARTIAL PIVOTING: Меняем R1 и R2 местами (pivot = 1 >> 0.001): 1·x + 1·y = 2 0.001·x + 1·y = 1 Теперь множитель = 0.001/1 = 0.001 - числа маленькие, ошибка минимальна ВСЕ ПРОМЫШЛЕННЫЕ РЕАЛИЗАЦИИ (numpy, LAPACK, MATLAB) используют pivoting по умолчанию. Вручную реализовывать без него - ошибка.

**Число обусловленности** (condition number) - метрика того, насколько система чувствительна к погрешностям. `np.linalg.cond(A)` больше 10^6 - система плохо обусловлена, решение ненадёжно. В машинном обучении плохо обусловленные матрицы встречаются при мультиколлинеарности признаков. Решение: регуляризация (ridge) или pseudo-inverse через SVD.

Сложность и когда метод Гаусса не подходит

Сложность и когда метод Гаусса не подходит

Для больших разреженных матриц (sparse) - тысячи уравнений, большинство коэффициентов нуль - промышленность использует **итеративные солверы**: Conjugate Gradient, GMRES, BiCGSTAB. Они не строят треугольную форму явно, а приближаются к решению итерационно. Метод Гаусса - основа теории, LU-разложение - промышленный стандарт для плотных матриц.

Метод Гаусса в реальных системах

Прямые и косвенные наследники алгоритма

Практика: метод Гаусса

Практика: метод Гаусса

Вопросы для собеседования

Что унести из урока

• **Метод Гаусса = последовательное исключение**: прямой ход строит треугольную форму, обратный ход читает ответ
• **LU-разложение = Гаусс в матричном виде**: P·L·U, O(n³) один раз, потом O(n²) для любого b
• **np.linalg.solve под капотом** - именно LU с partial pivoting (LAPACK dgesv)
• **Нормальное уравнение регрессии** θ = (XᵀX)⁻¹Xᵀy решается через solve, не через inv
• **Partial pivoting** - обязателен для численной устойчивости, все промышленные реализации включают его
• **Cond(A) >> 1** - признак плохо обусловленной системы; решение: ridge regularization или lstsq через SVD
• **Гаусс-Жордан** расширяет метод до нахождения обратной матрицы: [A|I] -> [I|A⁻¹]

Куда дальше

Метод Гаусса - фундамент для более мощных инструментов

• Обратная матрица — Гаусс-Жордан - стандартный алгоритм вычисления обратной
• Определитель матрицы — det вычисляется за O(n³) через треугольную форму (= произведение диагональных элементов U)
• SVD и псевдообратная матрица — Более устойчивая альтернатива для плохо обусловленных систем и прямоугольных матриц
• Линейная регрессия — Нормальное уравнение - прямое применение LU-разложения

Связанные уроки

• nm-01
• nm-01