Некоторые линейные функции обратимы: например, поворот на угол. Другие — необратимы: например, проекция. Понятие обратимости можно продолжить и на матрицы.
Определение. Матрица $A$ является обратимой, если существует матрица $A^{-1}$ такая, что
$$ A \cdot A^{-1} = A^{-1} \cdot A = I $$Из свойств матричного умножения следует, что для того, чтобы обратная матрица существовала, матрица $A$ должна быть квадратной, но этого не всегда достаточно.
Системы линейных уравнений
Понятие обратимости матриц важно для решения линейных уравнений:
$$ Ax = b \leftrightarrow \begin{cases} a_{11} x_1 + a_{12} x_2 + \ldots + a_{1n} x_n = b_1 \\ a_{21} x_1 + a_{22} x_2 + \ldots + a_{2n} x_n = b_2 \\ \ldots \\ a_{n1} x_1 + a_{n2} x_2 + \ldots + a_{nn} x_n = b_n \end{cases} $$ Если матрица $A$ обратима, то решение будет единственным, а именно $$ x = A^{-1} b $$В противном случае будет либо ноль, либо бесконечное количество решений, в общем случае живущих на каком-то многомерном пространстве. Об этом удобно думать так: каждое условие (строчка-уравнение) задает какую-то гиперплоскость в $n$-мерном пространстве, и пересечения этих гиперплоскостей либо пустые, либо совпадают, либо дают какое-то подпространство меньшей размерности.
Решение систем линейных уравнений — одна из основных задач вычислительной линейной алгебры, часто появляющаяся как подзадача в других алгоритмах.
Метод Гаусса
Метод Гаусса (англ. Gaussian elimination) позволяет решать системы линейных уравнений, а также находить обратные к матрицам и решать смежные задачи.
Он основывается на том факте, что с системами уравнений мы можем свободно делать следующие операции:
- поменять два уравнения местами,
- домножить любое уравнение на ненулевую константу,
- вычесть из одного уравнения другое.
При этом все вышеприведенные операции обратимы и не теряют никакую информацию, потому что все строки всегда будут нетривиальными линейными комбинациями друг друга, и любое изначальное уравнение всегда будет учтено с каким-то коэффициентом в каком-нибудь из уравнений.
Цель алгоритма — используя эти три операции привести систему к виду
$$ \begin{cases} 1 \cdot x_1 + 0 \cdot x_2 + \ldots + 0 \cdot x_n = b_1' \\ 0 \cdot x_1 + 1 \cdot x_2 + \ldots + 0 \cdot x_n = b_2' \\ \ldots \\ 0 \cdot x_1 + 0 \cdot x_2 + \ldots + 1 \cdot x_n = b_n' \end{cases} $$и тогда решением будет просто $x = b'$.
Алгоритм
Будем действовать в предположении, что решение существует и единственно.
Припишем вектор $b$ справа от матрицы $A$, и будем итеративно приводить получившуюся $n \times (n+1)$ матрицу к требуемому виду. На шаге $i$, мы хотим сделать так, чтобы в ячейке $a_{ii}$ стояла единица, а во всех остальных ячейках $i$-того столбца стояли нули.
Чтобы выполнить $i$-тый шаг:
- найдем какой-нибудь ненулевой элемент на $i$-том столбце и поменяем его строку местами с $i$-той;
- поделим всю $i$-тую строку на элемент $a_{ii}$ (он ненулевой: мы его специально искали на предыдущем шаге);
- пройдемся по всем остальным строкам $j \ne i$ и вычтем $i$-тую строку, помноженную на коэффициент $a_{ji}$.
В результате после $i$-того шага элемент $a_{ii}$ будет единичным (потому что мы его поделили на самого себя), а во всех остальных ячейках $i$-того столбца стали нули (потому что мы в этих позициях вычли $a_{ij} \cdot a_{ii} = a_{ij}$).
Таким образом, через $n$ шагов мы приведем матрицу к нужному виду, и ответом будет последний, $(n+1)$-ый приписанный столбец.
Реализация
Из соображений вычислительной стабильности имеет смысл на каждом шаге брать не любой ненулевой элемент в качестве опорного, а самый большой — мы ведь на него делим, а в случае в числами с плавающей точкой могло получиться очень близкий к нулю число.
typedef vector< vector<double> > matrix;
vector<double> gauss(matrix a) {
int n = (int) a.size();
for (int i = 0; i < n; i++) {
// находим опорный элемент
int best = i;
for (int j = i + 1; j < n; j++)
if (abs(a[j][i]) > abs(a[best][i]))
best = j;
// меняем строчки местами
swap(a[best], a[i]);
// нормализуем i-тую строчку (слева нули, с ними ничего делать не надо)
for (int j = n; j >= i; j--)
a[i][j] /= a[i][i];
// зануляем все остальные ячейки в i-том столбце
for (int j = 0; j < n; j++)
if (j != i)
// слева только нули, так что можно пройтись только по правой части
for (int k = n; k >= i; k--)
a[j][k] -= a[i][k] * a[j][i];
}
vector<double> ans(n);
for (int i = 0; i < n; i++)
ans[i] = a[i][n];
return ans;
}
Алгоритм работает $O(n^3)$.
Также его можно обобщить не только на действительные числа, но и на другие поля — например, на остатки по простому модулю, только нужно использовать соответствующую процедуру для деления.
Булевы матрицы
В контексте олимпиад, в большинстве задач, где требуется решать системы линейных уравнений, на самом деле это нужно делать над полем $\mathbb{Z}_2$, то есть по модулю 2.
Задача. Есть $n$ лампочек и $n$ переключателей: каждый активированный переключатель меняет состояние (включает или выключает) какого-то подмножества лампочек. Известно состояние всех лампочек, нужно восстановить состояние переключателей.
Нас по сути просят решить следующую систему:
$$ \begin{cases} a_{11} x_1 + a_{12} x_2 + \ldots + a_{1n} x_n \equiv b_1 \pmod 2\\ a_{21} x_1 + a_{22} x_2 + \ldots + a_{2n} x_n \equiv b_2 \pmod 2\\ \ldots \\ a_{n1} x_1 + a_{n2} x_2 + \ldots + a_{nn} x_n \equiv b_n \pmod 2 \end{cases} $$Здесь $x$ — состояния переключателей, $b$ — состояния лампочек, $A$ — информация о том, влияет ли переключатель на лампочку.
В таком случае можно значительно ускорить и упростить обычный метод Гаусса через битсет:
typedef bitset<maxn> t;
typedef array<t, maxn> matrix;
t gauss(matrix a) {
for (int i = 0; i < n; i++) {
int nonzero = i;
for (int j = i + 1; j < n; j++)
if (a[j][i])
nonzero = j;
swap(a[nonzero], a[i]);
for (int j = 0; j < n; j++)
if (j != i && a[j][i])
a[j] ^= a[i];
}
t x;
for (int i = 0; i < n; i++)
x[i] = a[i][n] ^ a[i][i];
return x;
}
Обратим внимание, что как в этой, так и в предыдущей реализации предполагалось, что решение существует и единственное. Если это не так, то в какой-то момент мы не найдем ненулевой элемент в качестве опорного.
В этом случае можно дальше запустить алгоритм, игнорируя полностью нулевые столбцы, и тогда в конце будет какое-то количество нулевых сумм, равных каким-то $b_i'$. Если все эти $b_i'$ нулевые, то решение существует, в противном случае — нет.
В случае входных данных из «реального мира» можно добавить небольшой шум $\epsilon_{ij}$ к коэффициентам: это гарантирует, что решение, причем уникальное, всегда будет существовать.
Базисы
Определение. Базисом множества векторов называется набор векторов, через линейную комбинацию которых единственным способом можно выразить все вектора этого множества и только их.
Базисы есть не только в линейной алгебре. Например, ${1, x, x^2}$ является базисом всех квадратных трёхчленов. Или ${\neg, \land, \lor}$ является базисом всех логических выражений (то есть всё можно выразить через И, ИЛИ и НЕ). В произвольном языке программирования можно выделить какой-то набор команд, через который можно будет написать всё что угодно, и он тоже в этом смысле будет базисом.
Определение. Система векторов ${a_1, a_2, \dots, a_n}$ называется линейно зависимой, если существует такой набор действительных коэффициентов $(\lambda_1, \lambda_2, \dots, \lambda_n)$, что $\lambda_i \neq 0$ хотя бы для одного $i$ и
$$ \lambda_1 \cdot a_1 + \lambda_2 \cdot a_2 + \dots + \lambda_n \cdot a_n = 0 $$В противном случае система называется линейно независимой.
Утверждение. Базис $n$-мерного пространства — это линейно независимый набор из $n$ векторов.
Доказательство. Пусть базис линейно зависим. Тогда вектор нулевой длины можно выразить хотя бы двумя способами, что противоречит определению.
Базисы часто требуется искать или поддерживать в олимпиадных задачах — часто вида «набрать базис максимального веса, где каждый вектор сколько-то стоит». Для проверки, является ли набор векторов базисом, можно проверить, выражается ли как-то нетривиально с помощью них нулевой вектор — это можно сделать методом Гаусса.
Упражнение. Дан массив из $10^5$ целых чисел от $0$ до $(2^{30}-1)$. Требуется найти количество различных подпоследовательностей этого массива, xor-сумма которых равна заданному числу $x$.