Матрицы и рекуррентности

Между 21 аэропортом летают рейсы: из Алматы в Астану их выполняют 72 разных позывных, из Алматы в Балхаш — один. Сколько разных цепочек ровно из k рейсов ведут из одного аэропорта в другой, если каждый позывной считать отдельным способом? При k = 1 ответ прочитан в таблице рейсов. При k = 2 его надо собрать по всем промежуточным аэропортам. А при k = 10 цепочек больше, чем помещается в 64-битное число. Нажмите на аэропорт и двигайте k.

Все эти числа — строка одной матрицы: таблицы рейсов C, возведённой в степень k. Из Алматы в Астану двумя рейсами — 1192 цепочки, а обратно в Алматы — 7358: улететь и вернуться можно через любой город, куда из Алматы летают. Степень считается быстрым возведением: C¹⁰ = C²·C⁸, и вместе с возведениями в квадрат это четыре произведения матриц вместо девяти. Переключитесь на «куда занесёт»: путешественник каждый раз улетает случайным рейсом. После первого перелёта он точно не в Алматы, после второго — с вероятностью 51 % снова там: почти все рейсы связывают города с двумя главными аэропортами, и путешественника качает между ними. Но качка затухает, и через десяток перелётов уже неважно, откуда он вылетел.

Идея

Пусть C_ij — число способов сделать один шаг из i в j. Тогда способов дойти за два шага — сумма по промежуточной точке t произведений C_it·C_tj, а это ровно элемент произведения матриц C·C. За k шагов — элемент C^k. Матрица шага бывает таблицей рейсов, таблицей переходов автомата или правилом, по которому из прошлых членов последовательности получается следующий; степень за O(log k) произведений даёт быстрое возведение в степень — ему всё равно, умножать числа или матрицы.

Матрицы смежности и произведения матриц — это главы о графах и о линейных отображениях «Царицы наук», а степени через собственные векторы, формула Бине для чисел Фибоначчи и цепь Маркова с самокатами — глава о собственных векторах. Здесь — как этим считать: какие задачи прячут в себе степень матрицы, сколько она стоит и когда её лучше не возводить.

Рекуррентность — тоже степень матрицы

Сколькими способами можно выложить полосу 3 × n доминошками 2 × 1? Идём по полосе слева направо, столбец за столбцом. Всё, что надо помнить о уже уложенной части, — какие клетки следующего столбца уже заняты лежачими доминошками, которые торчат из предыдущего. Это маска из трёх битов, восемь состояний. Матрица переноса T говорит, сколькими способами заполнить столбец, войдя в него в состоянии s и оставив после себя состояние t. Замощение полосы — путь из n шагов по состояниям от ровного края к ровному краю, и по утверждению выше таких путей (Tⁿ)₀₀.

2 x 10: 89 tilings, the matrix is 4 x 4
3 x 10: 571 tilings, the matrix is 8 x 8
4 x 10: 18061 tilings, the matrix is 16 x 16
8 x 8 : 12988816 tilings, the matrix is 256 x 256
4 x 10^18: 328515267 mod 1e9+7, 82 matrix products instead of 10^18 steps

Полоса шириной 2 замощается числом Фибоначчи способов: 89 для длины 10 — та самая рекуррентность F_n = F_n−1 + F_n−2, которую обычно и записывают матрицей 2 × 2. Здесь её никто не выписывал: матрица переноса построилась сама, перебором заполнений одного столбца. Шахматная доска 8 × 8 замощается 12 988 816 способами — это число в 1961 году нашли формулой Питер Кастелейн и, независимо, Харолд Темперли с Майклом Фишером; программа получает его перемножением матриц 256 × 256. Для длины 10¹⁸ шаги по одному заняли бы тысячи лет даже на быстром компьютере, а возведение в степень делает 82 произведения матриц по модулю.

Любая линейная рекуррентность порядка d — a_n = c₁a_n−1 + … + c_da_n−d — устроена так же: состояние — вектор из d последних членов, шаг — умножение на матрицу d × d, у которой в первой строке коэффициенты, а ниже единицы сдвигают вектор. Цена — O(d³ log n). Для рекуррентностей порядка в сотни и тысячи есть способ быстрее, через остатки от деления многочленов, — O(d² log n), а с быстрым умножением многочленов и того меньше; а если рекуррентность неизвестна, но известны первые члены, её восстанавливает алгоритм Берлекэмпа — Мэсси из статьи о генераторах случайных чисел.

Слова, которых нет: автомат вместо графа

Фермент рестрикции EcoRI режет ДНК там, где встречает шесть букв GAATTC. Сколько таких мест ждать в последовательности длиной 48 502 буквы, если буквы случайны? Среднее посчитать легко — каждая из 48 497 позиций начинает участок с вероятностью около 1/4096. Но чтобы узнать вероятность того, что участков не больше пяти, нужно сосчитать, сколько последовательностей их избегают, — а это пути в графе автомата Кнута — Морриса — Пратта. Его состояние — сколько букв участка уже совпало; каждая буква переводит автомат в новое состояние, и переходы задаёт та же префикс-функция, что и в поиске подстроки. Добавим в состояние счётчик найденных участков, взвесим переходы вероятностями букв, и строка матрицы в степени 48 502 — это распределение числа участков.

Геном фага лямбда — вируса, который заражает кишечную палочку, — содержит GAATTC всего пять раз, а случайная последовательность той же длины с теми же частотами букв — в среднем 11,9 раза; пять или меньше у неё выпадает с вероятностью около 2 %. С участком BamHI, GGATCC, картина та же: пять против 11,7. У митохондрии человека разница в пределах случайности, у маленького фага φX174 участков нет, но и ждали бы их всего полтора. Нехватку участков узнавания ферментов в геномах фагов заметили давно — в 1986 году об этом писал Пол Шарп, объясняя её отбором: у бактерий такие ферменты режут чужую ДНК. Но правило это не всеобщее, а модель «каждая буква независима» груба: в настоящих геномах соседние буквы связаны. Программа считает то же для трёх геномов и двух ферментов.

lambda phage         48502  EcoRI GAATTC: found 5, expected 11.86, P(none) = 6.96e-06, P(<= 5) = 0.0221
lambda phage         48502  BamHI GGATCC: found 5, expected 11.68, P(none) = 8.29e-06, P(<= 5) = 0.0246

Куда занесёт: цепь Маркова

Поделим каждую строку таблицы рейсов на её сумму: P_ij — вероятность, что путешественник из i улетит в j, если он берёт любой рейс аэропорта с равными шансами. Это цепь Маркова: куда лететь дальше, зависит только от того, где он сейчас. Распределение после k перелётов — строка P^k, и считается оно той же степенью или, если нужна одна стартовая точка, k умножениями вектора на матрицу. У нашей сети у этих строк есть предел: Алматы — 32,2 %, Астана — 22,8 %, независимо от того, откуда стартовать. И он почти совпадает с долей рейсов, которые вылетают из аэропорта: 31,8 % и 22,1 %. Это не совпадение.

Настоящая сеть почти симметрична — в одну сторону рейсов иногда на один-два больше, — поэтому и предел почти равен долям. Важнее, как быстро к нему подходят. Расстояние до предела (полусумма модулей разностей вероятностей) после первого перелёта 0,33, после второго 0,19, третьего 0,11 — каждый раз примерно в 0,57 раза меньше, и при этом вероятность Алматы прыгает: 0, 51 %, 22 %, 38 %, 29 %. Множитель 0,568 (программа ниже печатает его точнее) — модуль второго по величине собственного значения P, а то, что само это значение отрицательно, и даёт качку: сеть похожа на звезду, где из большого города летят в малые, а из малых — обратно в большие. Про собственные значения и степенной метод — глава о собственных векторах и статья о PageRank, где тот же предел цепи Маркова расставляет по важности страницы.

Код: цепочки и вероятности

 k  to Astana                       back to Almaty
 1  72                            0
 2  1192                          7358
 3  599774                        181406
 4  23398952                      61553540
10  38516548485362965772          49395526650223704252
24  a number of 48 digits      48 digits
from k = 10 some counts no longer fit in a signed 64-bit integer

Число цепочек растёт как степень самого большого собственного значения C — почти в сто раз за каждый рейс, — и уже при десяти рейсах не помещается в 64 бита. Python считает его длинными числами сам, в Go для этого нужен math/big, а в олимпиадных задачах ответ просят по модулю: тогда каждое сложение и умножение берут по модулю, и числа остаются маленькими. Аккуратнее всего с произведением: два остатка до 10⁹ дают до 10¹⁸, это ещё помещается в 64 бита со знаком, а сумма нескольких таких произведений — уже нет, поэтому остаток берут после каждого слагаемого.

Вторая половина программы — случайный путешественник: вектор вероятностей умножается на P перелёт за перелётом, а предел получен двумя тысячами таких умножений.

12  0.323   0.227   0.00066
each flight multiplies the distance by about 0.568
the limit: Almaty 0.322, Astana 0.228; their shares of all callsigns: 0.318, 0.221

Сколько это стоит

Произведение двух матриц m × m — m³ умножений, степень k — не больше 2 log₂k произведений: O(m³ log k). Но часто нужна не вся матрица, а одна строка: цепочки из одного аэропорта, замощения от ровного края. Тогда можно идти шагами — вектор на матрицу, m² умножений за шаг, k·m² всего. Шаги выгоднее, пока k меньше примерно 2m·log₂k. Замер на сети из 21 аэропорта, по модулю 10⁹ + 7:

Шаги растут прямо пропорционально k, возведение — как его логарифм. На Python линии пересекаются между k = 64 и 256, на Go — около 256. При k = 10¹⁸ выбора нет — только степень. А если матрица большая, а k скромное, как в задачах на графы из тысяч вершин, возведение в степень проигрывает и по времени, и по памяти: m³ при m = 10⁴ — это 10¹² умножений на одно произведение.

Ещё одна цена — размер чисел. Без модуля элементы C^k растут на постоянное число цифр за шаг, и умножать их приходится длинной арифметикой; с модулем каждое число занимает одно машинное слово, но ответ — только остаток.

Где это работает

• Олимпиадные задачи. «Посчитайте по модулю 10⁹ + 7 при n до 10¹⁸» почти всегда значит линейную рекуррентность, пути в маленьком графе или слова без запретных кусков, то есть степень матрицы. Обе задачи ниже — из этого семейства.
• Цепи Маркова. Модели очередей, погоды, поведения покупателей, вероятностные модели текста; их распределение через k шагов — строка P^k, а предел ищут степенным методом, как в PageRank.
• Комбинаторика и физика. Метод матрицы переноса считает замощения, раскраски полос и разбиения, а в статистической физике — суммы по всем состояниям цепочек атомов.
• Биоинформатика. Сколько раз мотив должен встретиться в случайной последовательности и насколько неожиданно то, что он встретился реже, — та же матрица автомата в степени длины последовательности.

История

В 1913 году Андрей Марков разложил первые 20 000 букв «Евгения Онегина» на гласные и согласные и посчитал, как часто за гласной следует гласная, а за согласной согласная. Вышло, что буквы зависят от предыдущих, а доля гласных всё равно устойчиво держится около 43 %: закон больших чисел работает и для зависимых событий, если зависимость устроена как цепь. Так появился первый разобранный на настоящих данных пример цепи, которая теперь носит его имя. Метод матрицы переноса пришёл из физики: им считают суммы по состояниям цепочек атомов, а в 1961 году два способа сосчитать замощения прямоугольника доминошками — Кастелейна и Темперли с Фишером — сделали задачу о шахматной доске классикой. Сам приём «возвести матрицу шага в степень» никому не приписывают: это быстрое возведение в степень, применённое к ассоциативному умножению матриц.

Задачи

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