Математика для программиста Матрицы, линейные отображения, определители и системы уравнений
0%

Матрицы, линейные отображения, определители и системы уравнений

34 мин 7072 слов
#математика#линейная алгебра#матрицы#определитель#numpy

Матрицы, линейные отображения, определители и системы уравнений

Матрица — самый недопонятый объект в программистской математике. Её почти всегда объясняют как «двумерный массив чисел», после чего правило умножения выглядит как произвольный ритуал: «строку на столбец, зачем — не спрашивайте». Отсюда все дальнейшие беды: определитель кажется формулой из ниоткуда, метод Гаусса — школьной механикой, а np.linalg.inv — правильным способом решить систему (это не так).

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

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

Зачем это программисту: где матрицы уже есть в вашем коде

  • Графика, игры, робототехника. Каждый кадр GPU перемножает миллионы матриц 4×4: модель → мир → камера → проекция. Матрица — это способ склеить цепочку преобразований в один объект.
  • ML. Полносвязный слой — это y = W·x + b. Батч из 1024 примеров — одно матричное умножение. Обучение сети — это, по сути, программа, которая 10^18 раз выполняет A @ B.
  • Аналитика и БД. Ковариационная матрица XᵀX, регрессия, рекомендательные системы (матричная факторизация), а в GraphBLAS даже обход графа в ширину — это умножение разреженной матрицы на вектор.
  • Криптография и надёжность хранения. Erasure-коды в Ceph, HDFS, Backblaze — матрицы над конечным полем GF(256); восстановление данных = решение системы линейных уравнений.
  • Компиляторы. Перестановка и скос вложенных циклов в полиэдральной модели — это унимодулярные матрицы, действующие на решётке итераций.

Во всех пяти случаях вопрос один и тот же: что эта матрица делает с пространством и обратимо ли это.

Линейное отображение: строгое определение

Определение. Пусть V и W — векторные пространства над одним полем F. Функция T: V → W называется линейной (линейным отображением, гомоморфизмом векторных пространств), если для всех u, v ∈ V и всех c ∈ F:

(1) аддитивность:   T(u + v) = T(u) + T(v)
(2) однородность:   T(c·u)   = c·T(u)

Эквивалентная компактная форма — сохранение линейных комбинаций:

T(a·u + b·v) = a·T(u) + b·T(v)   для любых a, b ∈ F

Из (2) при c = 0 немедленно следует T(0) = 0: линейное отображение обязано оставлять начало координат на месте. Это первый практический фильтр. Функция f(x) = x + 5 не линейна, хотя её график — прямая. Функция f(x) = x² не линейна. Функция «нормировать вектор» не линейна. ReLU не линейна (в этом весь смысл нелинейностей в нейросетях: без них стопка слоёв схлопнулась бы в одну матрицу — см. трек Нейронные сети).

Почему линейность — жёсткое ограничение, и почему это хорошо

Линейных функций R → R ровно столько же, сколько вещественных чисел: каждая имеет вид f(x) = k·x. Это чудовищно бедный класс по сравнению со всеми функциями. Но именно бедность даёт всё остальное:

  • конечное описание: любое линейное T: R^n → R^m задаётся m·n числами, а не бесконечной таблицей значений;
  • предсказуемая сложность: применить — O(mn), обратить — O(n³), и это верно всегда, без «зависит от данных»;
  • композиционность: композиция линейных линейна, поэтому цепочку преобразований можно свернуть заранее;
  • точная теория ошибок: для линейных задач мы умеем доказывать оценки устойчивости (число обусловленности), чего почти нигде больше нет.

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

Аффинное — не линейное: откуда берутся однородные координаты

Сдвиг x ↦ x + t не линеен (нарушает T(0) = 0), но встречается постоянно: перенос объекта в сцене, bias в нейросети, смещение в регрессии. Отображение вида x ↦ A·x + t называется аффинным.

Трюк, спасающий всю компьютерную графику: поднять размерность на единицу. Вложим R³ в R⁴ как (x, y, z) ↦ (x, y, z, 1). Тогда аффинное преобразование становится линейным в R⁴:

[ A  t ] [ x ]   [ A·x + t ]
[ 0  1 ] [ 1 ] = [    1    ]

Это и есть однородные координаты. Практический выигрыш огромен: теперь поворот, масштабирование, сдвиг и даже перспективная проекция — объекты одного типа, и цепочка преобразований склеивается одним произведением матриц 4×4, которое GPU умеет делать в железе. Аналогичный трюк в ML: приписать к вектору признаков константу 1, и bias станет обычным столбцом весов.

Главная теорема: матрица — это образы базисных векторов

Теорема. Пусть B = (b₁, …, b_n) — базис V. Тогда для любого набора векторов w₁, …, w_n ∈ W существует ровно одно линейное отображение T: V → W с T(b_j) = w_j.

Доказательство (конструктивное). Любой v ∈ V единственным образом раскладывается как v = x₁b₁ + … + x_nb_n (единственность разложения — это определение базиса). Линейность вынуждает

T(v) = T(x₁b₁ + … + x_nb_n) = x₁·T(b₁) + … + x_n·T(b_n) = x₁w₁ + … + x_nw_n

Правая часть определена корректно и задаёт линейное отображение; никакого выбора у нас не было — значит, оно единственно. ∎

Эта теорема — вся суть матриц. Знать линейное отображение = знать, куда уехали базисные векторы. Всё остальное достраивается по линейности.

Определение матрицы отображения. Зафиксируем базис B в V (размерность n) и базис C в W (размерность m). Матрицей T в этих базисах называется таблица A размера m × n, j-й столбец которой — координаты вектора T(b_j) в базисе C.

Для стандартных базисов это читается совсем просто: j-й столбец A равен A·e_j.

Матрица как линейное отображение: столбцы — образы базисных векторов

Разбор руками

Возьмём A = [[2, 1], [1, 2]]. Первый столбец (2,1) — это куда уехал e₁ = (1,0). Второй столбец (1,2) — куда уехал e₂ = (0,1). Куда уедет (3, -1)? По линейности:

A·(3, -1) = 3·A·e₁ + (-1)·A·e₂ = 3·(2,1) − 1·(1,2) = (6−1, 3−2) = (5, 1)

Никакого «строка на столбец» — просто линейная комбинация столбцов с коэффициентами из x. Это и есть столбцовая картина умножения:

A·x = x₁·(столбец 1) + x₂·(столбец 2) + … + x_n·(столбец n)

Из неё сразу видно центральное следствие: множество всех A·x — это линейная оболочка столбцов, то есть col(A). Система A·x = b разрешима тогда и только тогда, когда b лежит в этой оболочке.

Строчная картина — тот же результат, прочитанный по-другому: i-я координата A·x есть скалярное произведение i-й строки на x. Эта картина полезна для геометрии систем (каждое уравнение — гиперплоскость) и для реализации: она подсказывает обход по строкам, дружелюбный к кэшу в row-major раскладке.

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

import numpy as np

A = np.array([[2.0, 1.0],
              [1.0, 2.0]])
x = np.array([3.0, -1.0])

# столбцовая картина: линейная комбинация столбцов
col_view = x[0] * A[:, 0] + x[1] * A[:, 1]
# строчная картина: скалярные произведения строк на x
row_view = np.array([A[i] @ x for i in range(A.shape[0])])

print(col_view, row_view, A @ x)   # [5. 1.] [5. 1.] [5. 1.]

# j-й столбец — это буквально образ j-го базисного вектора
I = np.eye(2)
assert np.allclose(A @ I[:, 0], A[:, 0])
assert np.allclose(A @ I[:, 1], A[:, 1])

Матрица — это координаты, а не сама функция

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

Матрица так же относится к отображению, как строковое представление числа — к самому числу: 255, 0xFF и 0b11111111 — одно значение в трёх «базисах».

Умножение матриц — это композиция отображений

Пусть T: U → V имеет матрицу B (размер n × p), а S: V → W — матрицу A (размер m × n). Какова матрица композиции S ∘ T?

Считаем по определению: j-й столбец искомой матрицы — это (S∘T)(e_j) = S(T(e_j)) = S(B·e_j) = A·(B·e_j). То есть j-й столбец результата равен A, применённой к j-му столбцу B. Расписав это покоординатно, получаем

(A·B)[i][j] = Σ_k A[i][k] · B[k][j]

Вот и всё. Правило «строка на столбец» — не соглашение, а теорема. Из этого вывода мгновенно объясняется всё, что обычно заучивают:

  • Согласование размеров. (m×n)·(n×p) = (m×p), потому что выход T должен быть входом S. Число столбцов первой = число строк второй — это условие «функции стыкуются».
  • Ассоциативность (AB)C = A(BC). Композиция функций ассоциативна по определению, значит и произведение матриц ассоциативно. Доказывать через суммы не нужно.
  • Некоммутативность AB ≠ BA. Сначала надеть носки, потом ботинки — не то же самое, что наоборот.
  • I — единица. Тождественное отображение.
  • Отсутствие сокращения. Из AB = AC не следует B = C, если A необратима: необратимая функция теряет информацию.

Матрицы n × n образуют кольцо с единицей и некоммутативным умножением, а обратимые — группу GL(n, F) (подробнее в статье Абстрактная алгебра).

import numpy as np

R = np.array([[0.0, -1.0],    # поворот на 90° против часовой
              [1.0,  0.0]])
S = np.array([[1.0,  1.0],    # сдвиг (шир) вдоль оси x
              [0.0,  1.0]])

print(R @ S)   # [[ 0. -1.]  сначала сдвиг, потом поворот
               #  [ 1.  1.]]
print(S @ R)   # [[ 1. -1.]  сначала поворот, потом сдвиг
               #  [ 1.  0.]]

# оба сохраняют площадь, но это разные преобразования
print(np.linalg.det(R @ S), np.linalg.det(S @ R))   # 1.0 1.0

Читать произведение A @ B @ C следует справа налево: первым к вектору применяется C. Это источник половины багов в графическом коде и в цепочках трансформаций признаков.

Транспонирование — это сопряжённое отображение

Aᵀ[i][j] = A[j][i] — определение тривиальное, смысл — нет. Ключевое тождество: для стандартного скалярного произведения

⟨A·x, y⟩ = ⟨x, Aᵀ·y⟩   для всех x, y

То есть Aᵀ — единственное отображение, которое «переносит A на другую сторону скалярного произведения». Отсюда:

  • (AB)ᵀ = BᵀAᵀ — порядок переворачивается, потому что композиция сопряжённых идёт в обратном порядке;
  • в обратном распространении ошибки градиент, идущий назад через слой y = W·x, умножается ровно на Wᵀ. Это не «удобное совпадение», а определение сопряжённого оператора: forward умножает на W, backward — на Wᵀ;
  • AᵀA всегда симметрична и положительно полуопределена: xᵀ(AᵀA)x = ‖Ax‖² ≥ 0. На этом факте держатся нормальные уравнения МНК и вся ковариационная статистика.

Сколько это стоит: от n³ до BLAS

Наивное умножение двух матриц n × n — три вложенных цикла: n³ умножений и n³ сложений, то есть 2n³ флопов, память O(n²). Сложность O(n³) при O(n²) данных означает арифметическую интенсивность O(n): на каждый загруженный из памяти байт приходится много арифметики. Именно поэтому матричное умножение — идеальная нагрузка для процессора и GPU, и именно поэтому его удаётся разогнать в десятки раз без изменения асимптотики.

Порядок циклов и кэш

Три цикла i, j, k можно расставить шестью способами. Арифметика одна и та же, скорость различается в разы:

# ijk: внутренний цикл читает строку A (последовательно) и СТОЛБЕЦ B (шаг n) — промахи кэша
# ikj: внутренний цикл читает строку B и пишет строку C — всё последовательно в row-major
for i in range(n):
    for k in range(n):
        a_ik = A[i][k]
        for j in range(n):
            C[i][j] += a_ik * B[k][j]     # A[i][k] в регистре, B и C идут подряд

Реальные библиотеки идут дальше: блочное разбиение (tiling) под размеры L1/L2/L3, упаковка панелей в непрерывную память, микроядро на SIMD-регистрах, распараллеливание по потокам. Каноническое описание этой архитектуры — статья Goto и van de Geijn «Anatomy of High-Performance Matrix Multiplication» (https://dl.acm.org/doi/10.1145/1356052.1356053); её идеи лежат в основе OpenBLAS и BLIS.

Практический вывод для инженера: никогда не пишите матричное умножение сами. Разрыв между наивным тройным циклом на Python и вызовом BLAS через numpy — три-четыре порядка. Проверьте на своей машине:

import numpy as np, time

n = 256
A = np.random.rand(n, n)
B = np.random.rand(n, n)

t = time.perf_counter()
C1 = A @ B                                  # уходит в BLAS (OpenBLAS / MKL / Accelerate)
t_blas = time.perf_counter() - t

t = time.perf_counter()
C2 = np.zeros((n, n))                       # честные три цикла на чистом Python
for i in range(n):
    for k in range(n):
        a = A[i, k]
        for j in range(n):
            C2[i, j] += a * B[k, j]
t_py = time.perf_counter() - t

print(np.allclose(C1, C2), f"BLAS={t_blas:.4f}s  Python={t_py:.2f}s  x{t_py/t_blas:.0f}")

Полезно помнить и про раскладку: numpy по умолчанию хранит массивы в C-order (по строкам), а LAPACK — фортрановская библиотека и ждёт column-major. Иногда scipy.linalg тихо делает копию; флаг overwrite_a=True и np.asfortranarray позволяют этого избежать в горячем коде.

Быстрое умножение: Штрассен и показатель ω

В 1969 году Штрассен показал, что O(n³) — не предел: разбив матрицы на блоки 2×2, можно обойтись семью блочными умножениями вместо восьми, ценой лишних сложений. Рекурсия даёт O(n^log₂7) ≈ O(n^2.807).

Наименьшее известное значение показателя ω (такого, что умножение выполнимо за O(n^(ω+ε))) на 2024 год — примерно 2.3715 (Alman, Duan, Vassilevska Williams, Xu, Xu, Zhou, https://arxiv.org/abs/2404.16349). Нижняя граница — очевидное ω ≥ 2; где именно между 2 и 2.37 лежит истина, не знает никто.

Практический статус:

  • Галактические алгоритмы. Всё, что быстрее Штрассена, имеет константы, делающие их бесполезными при любых реальных n. Это чистая теория сложности (см. Теорию сложности).
  • Штрассен реален, но нишев. Даёт выигрыш примерно от n ≳ 1000, но заметно хуже по численной устойчивости: у него нет поэлементной оценки ошибки, только норменная. В LAPACK по умолчанию не используется.
  • Разреженность важнее асимптотики. Если в матрице nnz ненулей, умножение на вектор стоит O(nnz). Для графа социальной сети с n = 10⁹ и средней степенью 100 разница между O(nnz) и O(n²) — это разница между «работает» и «не существует».

Ранг и четыре фундаментальных подпространства

Определение. rank(A) — размерность пространства столбцов col(A). Фундаментальная теорема: размерность пространства строк равна размерности пространства столбцов, поэтому rank(A) = rank(Aᵀ) ≤ min(m, n).

Матрица m × n порождает ровно четыре подпространства, и их взаимное расположение — самая полезная схема во всей линейной алгебре (её популяризовал Гилберт Стрэнг):

Четыре фундаментальных подпространства матрицы

Теорема о ранге и дефекте (rank–nullity). Для A: R^n → R^m

dim ker(A) + dim im(A) = n        (дефект + ранг = размерность входа)

Интуиция бухгалтерская: на входе n степеней свободы. Часть из них отображение уничтожает (ker), остальные выживают и образуют образ. Ничего не теряется и не появляется.

Три эквивалентных прочтения в терминах кода:

  • Уравнения. A·x = b имеет решение ⟺ b ∈ col(A). Оно единственно ⟺ ker(A) = {0}. Иначе множество решений — аффинное многообразие x₀ + ker(A).
  • Данные. rank матрицы признаков — это число «настоящих» независимых признаков. Если rank < n, у вас мультиколлинеарность, и коэффициенты регрессии не определены однозначно.
  • Сжатие. Матрица ранга r хранится в r(m + n) числах вместо mn (как произведение U·V узких матриц). Это основа матричной факторизации в рекомендательных системах и LoRA-адаптеров в LLM.
import numpy as np
from scipy.linalg import null_space, orth

A = np.array([[1.0, 2.0, 3.0],
              [2.0, 4.0, 6.0],     # = 2 × первая строка
              [1.0, 1.0, 1.0]])

r = np.linalg.matrix_rank(A)                 # 2 — считается через SVD с разумным порогом
col = orth(A)                                # ортонормированный базис образа, m × r
ker = null_space(A)                          # ортонормированный базис ядра, n × (n − r)
left_ker = null_space(A.T)                   # левое ядро, m × (m − r)

print(r, col.shape, ker.shape, left_ker.shape)   # 2 (3, 2) (3, 1) (3, 1)

# rank–nullity: 2 + 1 = 3
assert r + ker.shape[1] == A.shape[1]
# ядро ортогонально строкам, левое ядро ортогонально столбцам
assert np.allclose(A @ ker, 0)
assert np.allclose(left_ker.T @ A, 0)

Обратите внимание на matrix_rank: он не ищет нули точно, а считает сингулярные числа и отбрасывает те, что меньше относительного порога. В арифметике с плавающей точкой «ранг» — это всегда решение с порогом, а не факт.

Определитель

Аксиоматическое определение

Определитель обычно вводят формулой, из-за чего он выглядит как заклинание. На самом деле его определяют свойствами, а формула — следствие.

Определение. det — единственная функция от n столбцов матрицы n × n, обладающая тремя свойствами:

(D1) полилинейность: линейна по каждому столбцу при фиксированных остальных
(D2) кососимметричность: перестановка двух столбцов меняет знак
(D3) нормировка: det(I) = 1

Из (D2) сразу следует: если два столбца равны, то det = 0 (перестановка ничего не меняет, но обязана менять знак ⟹ d = −d ⟹ d = 0 в поле характеристики ≠ 2). Отсюда — главное: det A = 0 ⟺ столбцы линейно зависимы ⟺ A необратима.

Геометрия: знаковый объём

(D1)–(D3) — это в точности аксиомы объёма. Значит:

|det A| — коэффициент, на который A умножает любой объём. Знак det A показывает, сохраняется ли ориентация.

Единичный квадрат площади 1 переходит в параллелограмм площади |det A| (см. SVG выше: det = 3, площадь выросла втрое). Куб — в параллелепипед объёма |det A|. Любая измеримая фигура — тоже: это ровно якобиан в формуле замены переменных в кратном интеграле.

Три немедленных следствия, которые больше не надо запоминать:

  • det(AB) = det(A)·det(B) — растянуть объём сначала в b раз, потом в a раз = растянуть в ab раз. Мультипликативность определителя, доказательство которой «в лоб» через суммы кошмарно, здесь очевидна.
  • det A = 0 ⟺ отображение схлопывает пространство в меньшую размерность (объём стал нулём) ⟺ информация потеряна ⟺ обратной функции нет.
  • det(A⁻¹) = 1 / det(A), det(cA) = cⁿ·det(A) (растянули все n осей).

Отрицательный определитель — это отражение: правая тройка векторов стала левой. В графике det < 0 у матрицы модели означает вывернутые наизнанку нормали и неправильный backface culling — типичный баг при отрицательном масштабе по одной оси.

Формула Лейбница и почему её нельзя считать

Из аксиом выводится явная формула:

det A = Σ_{σ ∈ S_n} sign(σ) · A[1][σ(1)] · A[2][σ(2)] · … · A[n][σ(n)]

Сумма по всем n! перестановкам, знак — чётность перестановки (комбинаторика перестановок разобрана в статье Дискретная математика и комбинаторика). Для n = 3 это шесть слагаемых — знакомое «правило Саррюса». Для n = 20 это 2.4·10^18 слагаемых: формула имеет теоретическую ценность и нулевую вычислительную.

Разложение по строке (формула Лапласа) не лучше: T(n) = n·T(n−1) — та же O(n!).

Разбор руками. Возьмём матрицу системы, которую будем решать ниже:

A = [  2   1  -1 ]
    [ -3  -1   2 ]
    [ -2   1   2 ]

Разложение по первой строке:

det A = 2·det[[-1, 2], [1, 2]] − 1·det[[-3, 2], [-2, 2]] + (−1)·det[[-3, -1], [-2, 1]]
      = 2·(−1·2 − 2·1) − 1·(−3·2 − 2·(−2)) − 1·(−3·1 − (−1)·(−2))
      = 2·(−4) − 1·(−2) − 1·(−5)
      = −8 + 2 + 5 = −1

Определитель ненулевой ⟹ система будет иметь единственное решение. Ниже мы получим то же −1 из LU-разложения за O(n³).

Как считают на практике

Реальный алгоритм: привести матрицу к треугольному виду методом Гаусса и перемножить диагональ.

det A = (−1)^(число перестановок строк) · U[0][0] · U[1][1] · … · U[n−1][n−1]

Работает, потому что элементарная операция «прибавить кратное одной строки к другой» не меняет определитель (следует из D1+D2), перестановка строк меняет знак, а определитель треугольной матрицы равен произведению диагонали (все остальные слагаемые Лейбница содержат нулевой множитель). Стоимость — те же O(n³), что и у LU-разложения, которое всё равно нужно для решения системы.

Практическая ловушка — переполнение. Определитель матрицы 500×500 со значениями порядка единицы легко оказывается 10^±300: он ведёт себя как произведение n чисел. Для матриц ковариации в статистике (логарифм правдоподобия гауссианы содержит log det Σ) считать log(det(A)) — прямой путь к -inf или nan. Правильный инструмент — slogdet, возвращающий знак и логарифм модуля отдельно:

import numpy as np

rng = np.random.default_rng(0)
M = rng.normal(size=(400, 400))

print(np.linalg.det(M))                  # часто inf или 0.0 — переполнение/потеря значимости
sign, logabsdet = np.linalg.slogdet(M)   # знак и log|det| по отдельности — устойчиво
print(sign, logabsdet)

# так считают log-правдоподобие многомерной гауссианы:
#   -0.5 * (x−μ)ᵀ Σ⁻¹ (x−μ) − 0.5·logdetΣ − (k/2)·log(2π)

Чего определитель НЕ говорит

Это место, где ошибаются чаще всего.

  1. det — плохой индикатор «почти вырожденности». det(0.01·I₁₀₀) = 10^-200, но матрица идеально обусловлена: делить на 0.01 абсолютно безопасно. Наоборот, у матрицы Гильберта 12×12 определитель тоже крошечный, и вот она действительно чудовищна. По величине det эти случаи неразличимы.
  2. Мера вырожденности — сингулярные числа и число обусловленности, а не определитель. Проверять if det(A) == 0 в коде с плавающей точкой — ошибка: почти никогда не сработает точно.
  3. det не нужен для решения систем. Правило Крамера (x_i = det(A_i)/det(A)) красиво и полезно теоретически (например, для доказательств про целочисленные решения), но требует n+1 определителей — O(n⁴) при устойчивом счёте — и численно хуже метода Гаусса. Используйте его только для n = 2 в формуле «на бумаге».
  4. det не аддитивен: det(A + B) ≠ det A + det B. Полилинейность — по столбцам, не по матрице целиком.

Обратная матрица и теорема об обратимости

A⁻¹ — матрица обратного отображения: A·A⁻¹ = A⁻¹·A = I. Для квадратной A размера n × n следующие утверждения эквивалентны (invertible matrix theorem — самый полезный список в курсе):

1.  A обратима
2.  det A ≠ 0
3.  rank A = n
4.  столбцы A линейно независимы
5.  строки A линейно независимы
6.  ker(A) = {0}, то есть A·x = 0 ⟹ x = 0
7.  A·x = b имеет решение для любого b
8.  A·x = b имеет единственное решение для любого b
9.  столбцы A образуют базис R^n
10. 0 не является собственным значением A
11. Aᵀ обратима
12. A представима как произведение элементарных матриц

Каждый пункт — отдельный способ проверить одно и то же свойство, и в разных задачах удобен разный. В доказательстве обычно берут (6) — его проще всего установить; в численном коде смотрят на (3) через SVD.

Почему обратную матрицу не считают

Стандартный совет численной линейной алгебры: если в коде встретился inv(A) @ b — это ошибка. Надо писать solve(A, b).

Причины:

  • Скорость. Явное обращение — примерно 2n³ флопов, плюс потом умножение 2n². LU-разложение с последующей подстановкой — 2n³/3, то есть втрое дешевле. При нескольких правых частях b фактор считается один раз, а каждая подстановка стоит O(n²).
  • Точность. Обращение накапливает ошибку и разрушает структуру. У ленточной, разреженной или треугольной матрицы обратная, как правило, плотная: обращение разреженной матрицы 10⁶×10⁶ физически не поместится в память, а решить с ней систему итерационным методом — рутина.
  • Устойчивость. Классическая оценка Хайема: вычисленное solve даёт малую невязку в смысле обратного анализа ошибок, а inv(A) @ b таких гарантий не имеет.
import numpy as np

rng = np.random.default_rng(42)
n = 800
A = rng.normal(size=(n, n)) + n * np.eye(n)   # хорошо обусловленная
b = rng.normal(size=n)

x1 = np.linalg.solve(A, b)          # LU с частичным выбором, 2n³/3
x2 = np.linalg.inv(A) @ b           # так делать не надо: ~3× дороже и менее точно

# сравниваем невязки, а не «правильность»: истинного x мы не знаем
print(np.linalg.norm(A @ x1 - b), np.linalg.norm(A @ x2 - b))

Единственные оправданные случаи явного обращения: матрица очень мала (2×2, 3×3 в графике), нужна сама обратная как объект (ковариация оценок в статистике — но и там честнее cho_solve), или вы пишете учебный пример.

Системы линейных уравнений

Геометрия множества решений

Для A·x = b возможны ровно три исхода, и никаких других:

Ситуация Условие Множество решений
Единственное решение b ∈ col(A), ker(A) = {0} точка
Бесконечно много b ∈ col(A), ker(A) ≠ {0} x₀ + ker(A) — аффинное подпространство размерности n − r
Нет решений b ∉ col(A) ∅ (в приложениях — ищем ближайшее по МНК)

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

Строчная картина даёт наглядную геометрию: каждое уравнение aᵢ·x = bᵢ — гиперплоскость в R^n, решение — их пересечение. Три плоскости в R³ могут пересечься в точке, по прямой, по плоскости или не пересечься вовсе (образовав «треугольную призму») — это ровно четыре сценария из таблицы.

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

Решаем систему:

 2x +  y −  z =   8
−3x −  y + 2z = −11
−2x +  y + 2z =  −3

Расширенная матрица и прямой ход (обнуляем поддиагональ):

[  2   1  -1 |   8 ]
[ -3  -1   2 | -11 ]      R2 += 1.5·R1
[ -2   1   2 |  -3 ]      R3 += 1.0·R1

[  2   1    -1  |  8 ]
[  0   0.5   0.5|  1 ]
[  0   2     1  |  5 ]    R3 -= 4·R2

[  2   1    -1  |  8 ]
[  0   0.5   0.5|  1 ]
[  0   0    -1  |  1 ]

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

−z = 1              ⟹ z = −1
0.5y + 0.5(−1) = 1  ⟹ y = 3
2x + 3 − (−1) = 8   ⟹ x = 2

Проверка подстановкой сходится. Заодно получен определитель: перестановок строк не было, диагональ треугольной матрицы 2 · 0.5 · (−1) = −1 — то самое −1, что мы посчитали разложением по строке.

Сложность. Прямой ход: ≈ 2n³/3 флопов, обратный ход: n². Память O(n²) (обычно in-place поверх исходной матрицы). Это фундаментальная стоимость решения плотной системы — и она не меняется с 1810-х годов.

Частичный выбор ведущего элемента: почему без него нельзя

Наивный Гаусс делит на диагональный элемент A[k][k]. Если тот равен нулю — деление на ноль. Если он просто маленький — катастрофа тише, но хуже: её никто не заметит.

[ 1e-20   1 | 1 ]
[ 1       1 | 2 ]

Точное решение ≈ x = (1, 1). Без пивотинга множитель равен 1/1e-20 = 1e20, вторая строка становится (0, 1 − 1e20 | 2 − 1e20). В double 1 − 1e20 и 2 − 1e20 округляются до одного и того же числа −1e20: единица и двойка теряются целиком. Получаем y = 1, затем x = (1 − 1)/1e-20 = 0. Ответ (0, 1) вместо (1, 1) — 100% ошибки, при том что матрица прекрасно обусловлена (cond ≈ 2.6).

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

import numpy as np

A = np.array([[1e-20, 1.0],
              [1.0,   1.0]])
b = np.array([1.0, 2.0])

print(np.linalg.solve(A, b))          # [1. 1.] — LAPACK делает пивотинг за вас
print(np.linalg.cond(A))              # ≈ 2.6 — матрица отличная, проблема была в алгоритме

Тонкость, о которой стоит знать: даже частичный выбор в худшем случае допускает рост элементов в 2^(n−1) раз (матрицы Уилкинсона). Полный выбор (по всей подматрице) устойчивее, но стоит O(n³) сравнений и почти не векторизуется, поэтому все библиотеки используют частичный: на практике патологии не встречаются, а цена нулевая. Подробный разбор — у Хайема, «Accuracy and Stability of Numerical Algorithms», гл. 9.

LU-разложение и переиспользование фактора

Метод Гаусса на самом деле вычисляет разложение P·A = L·U, где P — перестановка, L — нижняя треугольная с единицами на диагонали (множители), U — верхняя треугольная. Это не деталь реализации, а главный практический приём:

Разложить один раз:   2n³/3 флопов
Решить для каждого b:  2n²  флопов  (прямая + обратная подстановка)

Если у вас 1000 правых частей (типично для симуляций, фильтров, оптимизации), разница между «факторизовать один раз» и «звать solve 1000 раз» — это n/3 раз, то есть сотни. В scipy это lu_factor / lu_solve.

Своя реализация LU с пивотингом

Учебная реализация, полностью соответствующая тому, что делает LAPACK dgetrf (без блочности и оптимизаций):

import numpy as np

def lu_partial_pivot(A):
    """Находит P·A = L·U. Возвращает perm (P как массив индексов строк), L, U.
    Сложность: O(n³) по времени, O(n²) по памяти."""
    A = np.array(A, dtype=float)          # работаем на копии
    n = A.shape[0]
    perm = np.arange(n)
    L = np.eye(n)
    for k in range(n - 1):
        # частичный выбор: наибольший по модулю элемент в столбце k ниже диагонали
        p = k + int(np.argmax(np.abs(A[k:, k])))
        if abs(A[p, k]) < 1e-14:
            raise ValueError("матрица вырождена в пределах машинной точности")
        if p != k:
            A[[k, p], :] = A[[p, k], :]
            L[[k, p], :k] = L[[p, k], :k]      # уже вычисленные множители едут со строкой
            perm[[k, p]] = perm[[p, k]]
        for i in range(k + 1, n):
            m = A[i, k] / A[k, k]              # множитель, по модулю ≤ 1 благодаря выбору
            L[i, k] = m
            A[i, k:] -= m * A[k, k:]
    return perm, L, np.triu(A)

def lu_solve(perm, L, U, b):
    """Решает L·U·x = b[perm]. Сложность O(n²)."""
    y = np.zeros_like(b, dtype=float)
    bp = np.asarray(b, dtype=float)[perm]
    for i in range(len(bp)):                       # прямая подстановка (L с единицами)
        y[i] = bp[i] - L[i, :i] @ y[:i]
    x = np.zeros_like(y)
    for i in reversed(range(len(y))):              # обратная подстановка
        x[i] = (y[i] - U[i, i + 1:] @ x[i + 1:]) / U[i, i]
    return x

def perm_sign(perm):
    """Знак перестановки = (−1)^(число инверсий)."""
    p, s = list(perm), 1
    for i in range(len(p)):
        for j in range(i + 1, len(p)):
            if p[i] > p[j]:
                s = -s
    return s

A = np.array([[ 2.0,  1.0, -1.0],
              [-3.0, -1.0,  2.0],
              [-2.0,  1.0,  2.0]])
b = np.array([8.0, -11.0, -3.0])

perm, L, U = lu_partial_pivot(A)
assert np.allclose(L @ U, A[perm])                 # P·A = L·U
x = lu_solve(perm, L, U, b)
print(x)                                           # [ 2.  3. -1.]
assert np.allclose(x, np.linalg.solve(A, b))

# определитель бесплатно: det(P)·det(A) = det(L)·det(U) = произведение диагонали U
det = perm_sign(perm) * np.prod(np.diag(U))
print(det, np.linalg.det(A))                       # ≈ -1.0  -1.0

Обратите внимание, что perm_sign здесь квадратичный — для наглядности; в реальном коде знак просто копят как ±1 при каждой перестановке.

Какой солвер выбрать

Плотный LU — универсальный ответ, но далеко не всегда лучший. Дерево решений:

Обусловленность: когда «правильный» ответ бесполезен

Число обусловленности измеряет, во сколько раз задача усиливает ошибки входа:

cond(A) = ‖A‖ · ‖A⁻¹‖     (для спектральной нормы = σ_max / σ_min)

Основная оценка: если правая часть известна с относительной погрешностью δ, то относительная погрешность решения может достигать cond(A)·δ. Отсюда правило большого пальца:

При решении системы теряется примерно log₁₀(cond(A)) верных десятичных знаков. В double их всего ~16.

cond(A) = 10¹⁰ означает, что от решения останется шесть значащих цифр — и это не вина алгоритма. Хороший солвер даёт малую невязку ‖Ax − b‖, но малая невязка не означает малой ошибки: при большом cond вектор может быть далёк от истинного, идеально удовлетворяя уравнению почти-точно.

Матрица Гильберта: канонический кошмар

H[i][j] = 1/(i + j + 1) — безобидная на вид, чудовищно обусловленная: cond(H_n) растёт примерно как e^(3.5n).

import numpy as np
from scipy.linalg import hilbert

for n in (5, 10, 12, 14):
    H = hilbert(n)
    x_true = np.ones(n)
    b = H @ x_true
    x = np.linalg.solve(H, b)
    err = np.linalg.norm(x - x_true) / np.linalg.norm(x_true)
    res = np.linalg.norm(H @ x - b) / np.linalg.norm(b)
    print(f"n={n:2d}  cond={np.linalg.cond(H):.2e}  ошибка={err:.2e}  невязка={res:.2e}")

Картина, которую вы увидите: невязка остаётся на уровне 10⁻¹⁶ во всех случаях (алгоритм работает безупречно), а относительная ошибка при n = 14 доходит до порядка единицы — ответ полностью потерян. Это самый важный урок численной линейной алгебры: ‖Ax − b‖ ≈ 0 не является доказательством правильности.

Что с этим делать на практике:

  • всегда логировать np.linalg.cond для матриц, приходящих из данных;
  • масштабировать признаки: столбцы в разных единицах (рубли и доли) сами по себе раздувают cond на порядки, и простая нормировка часто чинит всё;
  • при больших cond переходить к регуляризации (ридж: (AᵀA + λI)x = Aᵀb) или к усечённому SVD — обе техники меняют задачу на близкую, но устойчивую;
  • не сравнивать матрицы с нулём: rank, matrix_rank, пороги на сингулярные числа.

Переопределённые системы и метод наименьших квадратов

В реальных задачах уравнений обычно больше, чем неизвестных: 10⁶ измерений и 20 параметров. Точного решения нет (b ∉ col(A)), и разумная замена — минимизировать ‖Ax − b‖².

Геометрия ответа читается прямо со схемы четырёх подпространств: нужно спроецировать b на col(A). Значит, невязка b − Ax должна быть ортогональна всем столбцам:

Aᵀ(b − Ax) = 0   ⟺   AᵀA·x = Aᵀb        (нормальные уравнения)

Невязка живёт в левом ядре ker(Aᵀ) — это его физический смысл.

Trade-off между двумя способами решать МНК:

Способ Стоимость Устойчивость Когда применять
Нормальные уравнения + Холецкий mn² + n³/3 cond(AᵀA) = cond(A)² — квадрат! m ≫ n, хорошо обусловленная задача
QR (Хаусхолдер) 2mn² − 2n³/3 cond(A) значение по умолчанию, так делает lstsq
SVD ~2mn² + 11n³ лучшая, даёт решение минимальной нормы ранг неполный или близок к неполному

Квадрат обусловленности у нормальных уравнений — не мелочь: задача с cond(A) = 10⁸ (терпимо) превращается в cond(AᵀA) = 10¹⁶ (полная потеря в double). Демонстрация на полиномиальной регрессии, где матрица Вандермонда быстро становится плохо обусловленной:

import numpy as np

rng = np.random.default_rng(0)
t = np.linspace(0, 1, 200)
y = np.exp(t) + 1e-9 * rng.normal(size=t.size)

deg = 14
A = np.vander(t, deg + 1)                       # матрица Вандермонда — печально известна

c_ne = np.linalg.solve(A.T @ A, A.T @ y)        # нормальные уравнения: cond в квадрате
c_qr, *_ = np.linalg.lstsq(A, y, rcond=None)    # QR/SVD внутри — устойчиво

print(f"cond(A)    = {np.linalg.cond(A):.2e}")
print(f"cond(AᵀA)  = {np.linalg.cond(A.T @ A):.2e}")
print(f"невязка NE = {np.linalg.norm(A @ c_ne - y):.3e}")
print(f"невязка QR = {np.linalg.norm(A @ c_qr - y):.3e}")

Невязка у нормальных уравнений будет заметно хуже, а при deg ≈ 20 они просто развалятся, тогда как lstsq продолжит работать. Мораль: scipy.linalg.lstsq / numpy.linalg.lstsq вместо ручного inv(A.T @ A) @ A.T @ y — не педантизм, а разница между рабочим и сломанным кодом.

Смена базиса и подобие

Если P — матрица перехода, чьи столбцы — новые базисные векторы, записанные в старых координатах, то одно и то же отображение имеет матрицы:

A  — в старом базисе
B = P⁻¹ · A · P  — в новом базисе

Читается справа налево: «перевести новые координаты в старые (P), применить отображение (A), вернуться в новые (P⁻¹)». Матрицы, связанные таким соотношением, называются подобными. Подобие — отношение эквивалентности, а классы эквивалентности — это и есть «отображения с точностью до выбора координат».

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

det(P⁻¹AP) = det(P⁻¹)·det(A)·det(P) = det(A)   — множители сокращаются

Это и есть постановка задачи следующей статьи: существует ли базис, в котором B = P⁻¹AP максимально проста — диагональна (собственные векторы), почти диагональна (жорданова форма) или диагональна в ортонормированном базисе (спектральная теорема, SVD).

Матрица как шесть разных объектов

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

Классы матриц со специальной структурой существуют ровно потому, что структура позволяет заменить универсальный O(n³) алгоритм на дешёвый и более устойчивый:

Где это в реальном коде

Графика и робототехника. Поза объекта — матрица 4×4 [[R, t], [0, 1]], где R ортогональна (RᵀR = I, det R = 1 для собственного поворота). Композиция звеньев манипулятора — произведение таких матриц (прямая кинематика). Практический нюанс: после сотен перемножений R перестаёт быть строго ортогональной из-за накопления ошибок, и её периодически реортогонализуют (через QR или SVD), иначе объекты начинают незаметно «раздуваться». Обратная к жёсткому преобразованию считается не через inv, а по формуле [[Rᵀ, −Rᵀt], [0, 1]] — на порядок дешевле и точнее.

Графы. Для матрицы смежности A элемент (A^k)[i][j] равен числу путей длины ровно k из i в j — прямое следствие определения произведения (суммируем по всем промежуточным вершинам). PageRank — это степенной метод для стохастической матрицы. Обход в ширину = умножение разреженной матрицы на вектор в булевом полукольце; на этом наблюдении построен стандарт GraphBLAS (https://graphblas.org/), позволяющий выразить алгоритмы из Теории графов как линейную алгебру над полукольцами. Замена (+, ×) на (min, +) превращает то же самое умножение в алгоритм Флойда–Уоршелла.

ML. Слой сети — Y = X·Wᵀ + b, где X — батч [batch, in]. Forward-проход — матричное умножение, backward — умножение на транспонированные матрицы. LoRA-дообучение LLM держится ровно на утверждении «поправка к весам имеет низкий ранг»: вместо ΔW размера d×d хранят B·A с узким внутренним измерением r, сокращая число параметров с d² до 2dr. Подробности — в треке Машинное обучение.

Коды и хранение. Линейный код над GF(2) задаётся порождающей матрицей G; кодирование — c = G·m, проверка — синдром H·c. Reed–Solomon в erasure-кодировании (Ceph, HDFS, Backblaze) использует матрицы над GF(256), у которых любая квадратная подматрица обратима (свойство Вандермонда/Коши). Именно это гарантирует, что любые k уцелевших шардов из n восстанавливают данные: восстановление — буквально решение системы k уравнений. Арифметика поля разобрана в Абстрактной алгебре.

Обратная сторона той же медали: линейность — это уязвимость. Шифр Хилла (c = K·m над Z₂₆) ломается известным открытым текстом за один шаг — достаточно собрать n пар и решить систему на K. Любой чисто линейный примитив взламывается линейной алгеброй, поэтому в AES между линейными слоями обязательно стоит нелинейный S-box.

Компиляторы. Преобразования вложенных циклов (перестановка, скос, реверс) описываются целочисленными матрицами, действующими на пространстве итераций. Легальны те, что унимодулярны — целочисленные с det = ±1, потому что у них обратная тоже целочисленная, а значит, преобразование биективно отображает решётку итераций на себя, не порождая дробных индексов. Это ядро полиэдральной модели в GCC Graphite и LLVM Polly (https://polly.llvm.org/).

Типичные заблуждения

  1. «Матрица — это таблица чисел». Таблица — представление. Объект — линейное отображение плюс выбор базисов. Отсюда все ошибки с транспонированием и порядком умножения.
  2. «Умножение матриц — странное правило». Это композиция функций, выведенная, а не постулированная. Если правило кажется произвольным, значит, вывод не был усвоен.
  3. «AB = BA, ну хотя бы иногда можно переставить». Нельзя. Перестановка допустима только при доказанной коммутируемости (например, диагональные в одном базисе, или A и A⁻¹).
  4. «det = 0 — надёжный тест на вырожденность». В float — практически бесполезный. Используйте matrix_rank или отношение σ_max/σ_min.
  5. «Маленький det = плохая обусловленность». Нет: det(0.01·I₁₀₀) = 10⁻²⁰⁰ при cond = 1. Это независимые величины.
  6. «Решаю систему через inv(A) @ b». Втрое дороже и менее точно. solve всегда.
  7. «Невязка мала, значит ответ верный». При большом cond — нет. Невязка мала почти всегда; ошибка может быть огромной.
  8. «Нормальные уравнения — нормальный способ сделать регрессию». Они возводят обусловленность в квадрат. Для чего-то серьёзнее учебного примера берите lstsq (QR/SVD).
  9. «Ранг — целое число, его можно точно посчитать». Только для точной арифметики. В float ранг — функция порога, и это фундаментально, а не недостаток библиотеки.
  10. «У произведения разреженных матриц разреженность сохранится». Часто нет: заполнение (fill-in) при LU-разложении разреженной матрицы может съесть всю память. Существуют целые алгоритмы переупорядочивания (AMD, nested dissection), борющиеся именно с этим.
  11. «След и определитель — примерно про одно». tr(AB) = tr(BA) всегда, а det(A+B) ≠ det A + det B никогда. У них разные алгебраические свойства и разные применения.
  12. «Ортогональная матрица — это про перпендикулярные строки». Про ортонормированные: единичной длины и попарно ортогональные. Матрица из ортогональных, но не нормированных столбцов не удовлетворяет QᵀQ = I.

Мини-итог

  • Линейное отображение сохраняет сложение и умножение на скаляр, а потому обязано оставлять ноль на месте; сдвиги — аффинны и линеаризуются подъёмом размерности (однородные координаты).
  • Линейное отображение однозначно определяется образами базисных векторов; матрица — это таблица координат этих образов, её j-й столбец равен A·e_j.
  • Умножение матриц — композиция отображений; отсюда следуют ассоциативность, некоммутативность и правило согласования размеров.
  • Транспонирование — сопряжение относительно скалярного произведения; поэтому backward-проход умножает на Wᵀ.
  • Ранг — размерность образа; rank + dim ker = n; четыре фундаментальных подпространства полностью описывают разрешимость и единственность.
  • Определитель — знаковый коэффициент изменения объёма, единственная полилинейная кососимметричная нормированная функция столбцов; det = 0 ⟺ необратимость. Но det не измеряет обусловленность и не нужен для решения систем.
  • Метод Гаусса = LU-разложение; частичный выбор ведущего элемента обязателен; фактор считается один раз и переиспользуется для многих правых частей.
  • cond(A) определяет, сколько верных цифр останется от ответа; малая невязка ничего не гарантирует.
  • В коде: solve вместо inv, lstsq вместо нормальных уравнений, slogdet вместо log(det(...)), matrix_rank вместо сравнения det с нулём.

Источники

Смежные темы трека: язык доказательств (эквивалентности, «тогда и только тогда») — в статье Математическая логика и доказательства; перестановки и их знак — в Дискретной математике и комбинаторике; поля, кольца и группа GL(n, F) — в Абстрактной алгебре; подробный разбор арифметики с плавающей точкой и накопления ошибок — в Численных методах.

Что дальше

Мы научились применять и обращать линейные отображения, но по-прежнему смотрим на них в случайно выбранном базисе. Следующий шаг — найти правильный базис: тот, в котором отображение выглядит как простое растяжение вдоль осей. Это собственные векторы и собственные значения, спектральная теорема, а для прямоугольных матриц — сингулярное разложение, из которого следуют PCA, псевдообратная, малоранговые приближения и честное определение числа обусловленности.

Читайте дальше: Собственные значения, SVD и матричные разложения.

Нашли неточность? Выделите фрагмент текста — рядом появится жучок.

Нужен разбор именно вашей ситуации?

Статья описывает общий случай. Если у вас частный — можно разобрать его отдельно, платно. А если не хватает целого материала, предложите тему: её оплачивают вскладчину, и она выходит открытой для всех.

Доска запросов